Logboek Transcriptomics
Jarno Jacob Duiker
2024-10-22
intro
Transcriptomics
Transcriptomics is de analyse van de RNA transcriptie’s geproduceerd door het genotypen op een bepaald moment wat een link geeft tussen het genoom, het proteoom en het cellulaire phentype. Het is een stap in de Omics pipeline waar dingen zoals ziekten en genetische afwijkingen en het effect er van word onderzocht. Transcriptomics in ons onderzoek word gedaan door de RNA transcripties te vergelijken van 2 groepen muizen. Een jonge groep muizen van 8-10 weken en een oudere groep muizen van 100-120. In ons onderzoek kijken naar de expressie in een paar gen pathways: 04662, 04660, 04213, 04211, 04640.
De 04662 is de B cell receptor signaling pathway deze laat zien ““.
De 04660 is de T cell receptor signaling pathway deze laat zien ““.
De 04211 is de Longevity regulating pathway deze laat zien ““.
De 04640 is de Hematopoietic cell lineage deze laat zien ““.
In de gene pathways kan er worden gekeken naar de verandering in expressie, hier is dus te zien of de gevonden variant mutatie effect heeft op de expressie van het gen. Een gen pathway ziet er uit als volgt.
Hier is een door ons met KEGG + pathview library gemaakte gene pathway te zien. Rechts boven is een kleuren schaal te zien met -1, 0 en 1. De kleuren schaal en cijfers betekenen ook wat, wanneer de gen een witte/grijze kleur heeft is de verandering in de expressie 0 wat dus betekent dat de expressie van dit gen normaal is. Wanneer een gen block een groene kleur heeft en dus richting de -1 gaat zegt dit dat de expressie is minder is dan de norm. Wanneer een gen een rode kleur heeft en dus richting de 1 gaat zegt dit dat het een overexpressie heeft meer dan het normaal.
De verschillen in deze expressies hebben veel invloed op de werking van de pathway.
Hier is een flowchart te zien die laat zien hoe wij dit project gaan doen. We beginnen dus met de data downloaden, controleren de kwaliteit en trimmen. Wanneer de kwaliteit goed is gaan we readmappen met HISAT2 daarna data met samtools sorteren en indexeren. De Data wordt dan verwerkt met Feature Counts om de gen expressie vast te stellen. Dan word de data geanalyseerd op 2 manieren. Met de wilcoxon rank sum test, deze test beoordeelt niet alleen het teken van de verschillen, maar bepaalt ook de grootte van de waargenomen verschillen in de gegevens. De andere test is DEseq2, Deseq gebruikt een statistisch model om uitterekenen wat de verschillen in gen expressie zijn tussen 2 of meer groupen samples (RNA). Aan de hand van de feature count data en de DEseq2 uitkomt worden visualisaties gemaakt.
KEGG - is een knowledge base voor systematische anlyse van gen functies, waar genomics informatie gelinked kan worden met functionele info. Er zijn ook visualisatie die hier uit kunnen komen zoals pathways waar over en onder expressie van genen word laten zien in een gen pathway
Volcano plot - is een plot gemaakt door het plotten van de negatieve logartime van de p value op de y axis meestal met een base van 10. Met de log fold change op x axis. Een volcano plot laat de statistische significantie zien met P value tegen de groote van de fold change.
Hier zijn alle library’s die wij gebruiken in dit project.
## Loading required package: S4Vectors## Loading required package: stats4## Loading required package: BiocGenerics##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff, table,
## tapply, union, unique, unsplit, which.max, which.min##
## Attaching package: 'S4Vectors'## The following object is masked from 'package:utils':
##
## findMatches## The following objects are masked from 'package:base':
##
## expand.grid, I, unname## Loading required package: IRanges## Loading required package: GenomicRanges## Loading required package: GenomeInfoDb## Loading required package: SummarizedExperiment## Loading required package: MatrixGenerics## Loading required package: matrixStats##
## Attaching package: 'MatrixGenerics'## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.##
## Attaching package: 'Biobase'## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedians## Loading required package: ggplot2## Loading required package: ggrepel##
## Attaching package: 'data.table'## The following object is masked from 'package:SummarizedExperiment':
##
## shift## The following object is masked from 'package:GenomicRanges':
##
## shift## The following object is masked from 'package:IRanges':
##
## shift## The following objects are masked from 'package:S4Vectors':
##
## first, second##
## Attaching package: 'tidyr'## The following object is masked from 'package:S4Vectors':
##
## expandCode - path melding
Wanneer er geen input path staat bij het BASH script ga er dan vanuit dat de code in de map is uitgevoerd waar de input bestanden staan.*
SRA download
De SRA’s van de RNA seq worden gedownload met het onderstaande stuk code. De gebruikte accention list staat in de github repository.
De eerste stap is het downloaden van de genetische data die verkregen was door het sequencen van de samples. Met behulp van de SRA run selector van NCBI is een selectie van SRA’s gemaakt die gedownload moeten worden.
Een SRA (Sequence Read Archive) is een gecomprimeerd archief dat de sequencing reads bevat.
Om deze bestanden te downloaden wordt gebruikt gemaakt van “prefetch”, deze commandline tool is onderdeel van de SRA toolkit.
In de volgende stap worden deze bestanden uitgepakt.
prefetch $(</students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt) \
--output-directory "/students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/" --max-size 200GGebruikte samples:
In ons onderzoek gebruiken we de 2 muizen groepen, de oude en jonge muizen. De oude muizen zijn 100-120 weken oud en de jonge muizen zijn 8-10 weken oud. Deze twee groepen hebben allebei vier replicaten. De samples van alle muizen zijn verkregen met een Illumina NextSeq 500.
Oude muizen: (Genotype: Rag2-/-)
SRR21754423
SRR21754422
SRR21754421
SRR21754420
Jonge muizen: (Genotype: Rag2-/-)
SRR21754408
SRR21754417
SRR21754418
SRR21754419
Kwaliteits check
De gedownloade bestanden zitten nog in hun SRA format en dit is niet de filetype die wij gebruiken voor onze tools.
Als eerste moet de kwaliteit van de data gecheckt worden en dit wordt gedaan met het programma FastQC dit programma leest de ASCII characters uit de file die bij elke sequentie staan waardoor het weet wat de kwaliteit is van de data. De voledige werking van FastQC word later nog wat dieper uitgelegt.
Om de bestanden van .SRA om te zetten naar .fastq word fasterq-dump gebruikt. Deze kijkt in de SRA files en pakt ze uit naar hun fastq format. Fasterq-dump is deel van de SRA toolkit (onze gebruikte versie is 3.1.1)
Een FastQ file ziet er als volgt uit.
De file heeft een title waar de sample naam word weergegeven + de lengte. Daarna is de sequentie te zien, Dan is de titel te zien van de bijbehoorende kwaliteit sequentie. De laatste regel is dus de kwaliteits score sequenties deze is in ASCII
De kwaliteit word dus gescoord in ASCII characters, in het figuur hier onder is de score te zien van een voorbeeld sequentie, ook is te zien dat de score van het ! teken de laagst mogelijk kwaliteit is. De K is de hoogste kwaliteit mogelijk.
De formule die hier wordt gebruikt is de Phred kwaliteits score, deze score is logaritmisch gerelateerd aan de waarschijnlijkheid dat de “base call” verkeerd is.
Q = -log(E) waarbij Q de phred score is en E de waarschijnlijkheid van verkeerde base call.
De illumina gesequencenste sequenties hebben een afwijking van 33 hier word dus altijd -33 gedaan voordat de score wordt berekent.
Bron: https://gencoded.com/index.php/2020/05/20/fastq-format-an-overview/
Bron: https://gatk.broadinstitute.org/hc/en-us/articles/360035531872-Phred-scaled-quality-scores
# Cat leest elke File-parameter in volgorde van Acc_list_transcriptomics.txt en daarmee kijkt fasterq-dump welke files het moet uitpakken. Met parallel worden meerder files tegelijk verwerkt
cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | \
parallel fasterq-dump -O /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/FASTQ {}
# hierboven wordt de map output neergezetKwaliteitscontrole met fastqc
Voor het trimmen: FastQC wordt gebruikt voor het trimmen. Dit programma kan de scores lezen van de fastQ files en geeft hier een visuele output van. Dit zijn verschillende grafieken en tabellen hier onder is een tabel met uitleg voor elke visualisatie / tabel.
| Catagorie | Uitslag | Hoe lees je dit |
|---|---|---|
| Basic s tatistics | Hier krijg je een tabel met verschillende dingen zoals de Filename, file type, sequence length, %GC, Total sequences, Encoding en de sequences flagged as poor quality | Het aflezen van deze tabel is lezen en checken of er geen fouten in zitten. |
| Per base sequence quality | Een grafiek zoals te zien is in figuur 1 | De grafiek heeft 3 vlakken die van verschillende grote en kleuren. De kwaliteit ranged van 0-40, 28-40 is goed and gekleurd groen, 20-28 is gekleurd geel/oranje, 0-20 gekleurd rood. Verder representeert de gele box de 25th to 75th percentile. De zwarte lijnen geven de 10th en 90th percentile weer. De blauwe lijn geeft de gemiddelde scores voor kwaliteits controle score voor de nucleotide. gebaseerd op deze dingen is te zien dat de base van 1 tot 100 goede kwaliteit hierna gaat de kwaliteit sterk naar beneden. wat al aan geeft dat er getrimmed moet worden. |
| Per tile sequence quality | Hier is een plot te zien waar de reads worden weergegeven als een soort heatmap. waarbij elke tegel 1 een read weergeeft op een bepaalde positie deze kwaliteit van de reads worden vergeleken met elkaar. | Het aflezen van deze grafiek is door te kijken naar de kleur per tegel per positie. een donkerblauwe tegel betekent goede kwaliteit en hoe lichter de tegel word hoe slechter de kwaliteit dus lichtblauw betekent een slechte kwaliteit. Op de Y-as staat dan in welke tile het is en de X-as welke positie. |
| Per sequence quality score | Een grafiek die je de gemiddelde kwaliteit op de x-as en de nummer van sequences met gemiddelde op de y-as. |
Het beste is wanneer de meeste reads een hoge gemiddelde kwaliteits score hebben en er geen grote dip in de grafiek is, dit betekent een lage kwaliteit. |
| Per base sequence content | Hier is een plot te zien met een y-as waar 0-100 aangegeven wordt en een x-as met de “positie in read (bp)” in de grafiek staan 4 lijnen met het percentage per base | Het volgende is uit de grafiek te halen er is significante variatie in de nucleotide distributie aan het bein van de reads positie 1-10. Dit zou kunnen zijn door de voorbereiding of de vooroordelen in het sequencen. A, T, C en G zijn niet gelijk gerepresenteerd. na de 10de positioe zijn de base wat gestabilizeerd wat aan geeft dat de sequence kwaliteit in de rest van de reads hoger zijn. Het afgekeurde kruis komt dus door de choas van 1-10. |
| Per sequence GC content | Per sequence GC content geeft weer een plot weer met 2 lijnen. Een blauwe lijn die de theoretische distrubutie aangeeft wat dus een richtlijn is, en de GC count per read wat dus de gelezen data is. | Het aflezen wordt door de twee lijnen vergelijken. het is ideaal als de twee lijnen overlappen of dichtbij elkaar liggen. Wanneer er meerdere pieken zijn die afwijken van de blauwe theoretische lijn kan dat betekenen dat er misschien sprake is van besmetting of sequencing errors. Er komt dus een rood kruis wanneer de GC abnormaal is vergeleken met de theoretische verwachting wat dus zegt da de algemene kwaliteit niet goed is. |
| Per base N content | De per base N content grafiek geeft de frequentie weer van “N” basecalls op elke positie in de reads. De “N” staat voor een onzekere of onbekende nucleotide, deze kon de sequencer niet identificeren als 1 van de base (A, T, C, G) | X-as geeft de positie weer in de reads, Y-as geeft het percentage van reads met een “N” base op elke positie. Een hoge waarde betekent dat de sequencer op die positie vaak onzeker was over welke nucleotide er aanwezig was. De verwachting is dat er een zeer laag percentage N’s in de sequence zit, de standaard hiervoor is <1%. afwijkende resultaten kan wijzen op problemen zoals slechte kwaliteit van de sequentie. Ze komen vaak voor aan het begin of eind van de reads. |
| Sequence length dis tribution | De sequence length distribution grafiek laat de verdeling van de lengtes van de reads zien. X-as geeft de lengte van de sequenties (in basenparen) en de Y-as toont aantal reads van die specifieke lengte. | Deze grafiek is belangrijk bij NGS omdat afwijkingen in de sequentie lengte kunnen wijzen op fragmentatie of technische fouten tijdens sequencering . Een ideaal resultaat is een scherpe piek op 1 specifieke lengte wanneer je 150 bp doelreadlengte hebt zou de meeste sequencing output op precies 150 bp moeten vallen wat een piek bij die lengte zou moeten opleveren. Als er meerdere pieken zijn of een vrede spreiding van lengtes kan dit beteken dat er sequencing fouten, slechte adaptertrimming of degradatie van het DNA-monster. |
| Sequence du plication levels | De Sequence duplication levels grafiek laat het percentage van sequences zien die meer dan 1x voorkomen. Duplicaties kunnen voorkomen door technische “artifacts” tijdens de sequencing en andere factoreren zoals P CR - amplificatie, en kunnen de diversititeit van de dataset verminderen. Dit heeft invloed op downstream analyses zoals genexpressie of va r iantdetectie. X-as laat het aantal keren zien dat de specifieke sequentie wordt
gedupliceerd Y-as laat het percentage van reads dat voorkomt met duplicatie niveau zien |
Het aflezen van de grafiek is door te kijken naar de twee lijnen in
de grafiek de Duplicated sequences (meestal een rode lijn) geeft het
percentage met unieke sequenties weer zonder correctie voor natuurlijke
duplicaties. Het geeft een beeld van hoeveel van de reads in de dataset
meerdere keren voorkomen zonder te differentiëren tussen technische en
biologische duplicates. in de grafiek wil je graag een scherpe daling
zien waarbij de meeste reads een duplicatie niveau van 1 hebben en het
percentage gedupliceerede sequenties daarna snel afneemt. De total sequences lijn (meestal een blauwe lijn) corrigeert voor verwachte natuurlijke duplicates deze laat zien hoe de duplicatie eruit zou zien zonder technische artefacten en biedt een eerlijker beeld van hoeveel sequenties overgedupliceerd zijn de de sequencering zelf. Als er een groot verschil is tussen de lijnen betekent dit dat er duplicatie is onstaan door technische factoren zoals PCR duplicatie in plaats van biologische oorzaak. Te hoge aantal dupolactie kan probleem zijn voor downstream analyses en je wil dus dat de lijnen dichtbij elkaar liggen. |
| Overre presented sequences | Dit is een tabel die sequenties van op zijn minst 20bp die vaker voorkomen dan 0.1% van de totale nummer van sequenties. In de tabel staat de sequenties u itgeschreven, de count, het percentage en de w a arschijnlijke bron | Wanneer er een abnormaliteit te zien was in de Per sequence GC content grafiek kan er in deze tabel worden gekeken om de bron te vinden. als het niet staat als een bekende adapter of “vector”, kan het helpen om de data te blasten o te identiteit te vinden in de tabel. |
| Adapter content | Een grafiek waar verschillende adapters staan | Het geeft aan of de sequenties adapterfragmenten bevatten en van verschillende apparaten. als deze aanwezig zijn is het te zien door af te lezen in de grafiek welke positie er zijn om ze vervolgens weg te trimmen. deze “adaptercontent” is er voor identificatie van de DNA |
Tabel 1: Uitleg fasta uitlezing
Het volgende commando is uitgevoerd in de FASTQ/ directory.
cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
fastqc -o /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/fastqc_output/voor_trimmen/ {}_1.fastq {}_2.fastq Vervolgens met met “multiqc .” in de output directory een multiqc rapport gemaakt van de resultaten van fastqc.
MultiQC wordt gebruikt om naar alle fastqc output files te kijken. Multiqc zoekt een directory voor FastQC output en compileert deze files zodat de algemene resultaten over meerdere files beter te begrijpen is. De grafieken worden bijvoorbeeld over elkaar gelegd om de verschillen te zien.
Na analyse van de fastQC is te zien dat onze sequence reads korte sequenties zijn van 150 bp. Deze corresponderen met een deel van een DNA sequentie. Het genoom kan niet in zijn geheel gelezen worden door technische beperking van de gebruikte illumina machines. Hoe langer de reads hoe meer kans dat er fouten gemaakt worden tijdens het sequencen.
daarom als er genoeg reads zijn die het laatste stuk bevatten kan er een stuk afgeknipt worden. Ook moeten de adapters verwijderd worden, adapters zijn korte stukken DNA die op de flow cell zitten als linkers. bron: https://www.lubio.ch/blog/ngs-adapters
Voor het trimmen wordt gebruik gemaakt van Trimmomatic ( versie 0.39 ) Trimmomatic gebruikt Maximum Information Quality Filtering om te bepalen of een read moet worden getrimmed + de paramenters van de gebruiker. Maximum Information Quality Filtering probeert de formule te bepalen bij welke informatie de hoogste kwaliteit heeft door verschillende eigenschappen van de informatie te combineren, zoals: Relevantie (afgeleid via de sigmoidfunctie), Correctie of betrouwbaarheid (door het product van correctieprobabiliteiten). Het doel is om de score te maximaliseren voor informatie die zowel relevant als zeer betrouwbaar is.
De door ons opgegeven parameters staan in het codeblock deze zijn.
MINLEN:40 dit zegt dat de minimale lengte van de reads 40 is.
SLIDINGWINDOW 4:20 het eerste nummer specificeert de grootte van de sliding winden en het tweede nummer is de gemiddelde read kwaliteit binnen de window van het eerste nummer. Een sliding window is dat er steeds in dit geval 4 basen bekeken worden en dat keer op keer.
ILLUMINACLIP is het path naar een bestand dat adapter sequenties bevat voor illumina adapters. De gebruikte adapter: TruSeq3-Se.fa:2:30:10 De 2:30:10 betekenen het volgende: 2 is de “seed mismatch”, het aantal mismatches dat toegestaan is in een sequentie die een adapter kan zijn.30 is de “palindrome clip threshold en 10 is de”simple clip threshold”, specificeert hoe accuraat de match tussen de adapter sequentie en de mogelijke adapter in de read.
Hier is een foto te zien van hoe Trimmomatic weet wanneer hij moet
knippen 
Trimmomatic neemt 3 scores, de length Threshold, coverage factor en de Error rate factor en vormt hier 1 totale score is en trimt hem op de piek zoals te zien in de grafiek hierboven.
Dit is gebaseerd op deze formule.
Deze formule berekent een score \(\text{Score}(l)\) door een combinatie van drie elementen: een sigmoid-achtige factor om de score op basis van \(l\) te dempen, een machtsverhouding van \(l\) met een parameter \(s\), en het product van correctiewaarden \(P_{\text{corr}}[i]\), die allemaal samen de uiteindelijke score beïnvloeden.
Bron: http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/TrimmomaticManual_V0.32.pdf
cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
'TrimmomaticPE -threads 80 ' \
'/students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/FASTQ/{}_1.fastq' \
'/students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/FASTQ/{}_2.fastq' \
'/students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/trimmomatic_output/paired/{}_forward_paired.fastq' \
'/students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/trimmomatic_output/unpaired/{}_forward_unpaired.fastq' \
'/students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/trimmomatic_output/paired/{}_rev_paired.fastq' \
'/students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/trimmomatic_output/unpaired/{}_rev_unpaired.fastq' \
'ILLUMINACLIP:/students/2024-2025/Thema05/3dconformatieChromatine/Trimmomatic/Trimmomatic-0.39/adapters/TruSeq3-SE.fa:2:30:10' \
'MINLEN:40 ' \
'SLIDINGWINDOW:4:20'Kwaliteitscontrole na trimming met fastQC
Na het trimmen van de reads kijken we weer naar de kwaliteit dit is om te kijken of we de parameters moeten aan passen van de trimmomatic of dat we met deze getrimmde data kunnen werken.
cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
fastqc -o /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/fastqc_output/na_trimmen/ {}_forward_paired.fastq {}_rev_paired.fastq Indexeren referentie genoom
Om te kunnen zien of onze reads verschillen met de normale expressie is het belangrijk dat wij een referentie genoom meenemen in onze analyse. De tools die wij gebruiken hebben tevens het referentie genoom nodig zodat ze functioneren.
Door middel van readmapping worden de reads dus vergeleken met het bekende genoom, in ons geval de mm39 het op dit moment nieuwste muis referentiegeoom. De reads worden er dus tegen aan gelegd zodat hun locatie in het genoom bekend wordt.
Het programma dat gebruikt gaat worden voor het alignen is HIsat2, deze heeft dus zoals eerder al gezegd een referentie genoom nodig maar deze moet dus geindexeerd zijn en deze moet eerst gemaakt worden.
Dit kan gedaan worden met Hisat2 build. Het referentie genoom mm39 word gebruikt ook wel met de bestand naam GCF_000001635.27, In het bestand met het referentie genoom zijn twee referentie genomen aanwezig GCF_000001635.27 en GCA_000001635.9/. Het verschil is dat GCF van Refseq is en GCA van GeneBank.https://www.ebi.ac.uk/training/online/courses/functional-genomics-ii-common-technologies-and-data-analysis-methods/rna-sequencing/performing-a-rna-seq-experiment/data-analysis/read-mapping-or-alignment/
Bron: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3375638/
Bron Dowload ref genoom : https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000001635.27/
In de onderstaande code chunk wordt het geïndexeerde muisgenoom “hisat2_index” genoemd.Daarna worden twee samples (forward read en reverse read) ge-aligned met het zojuist verkregen geïndexeerde referentiegenoom.Deze twee bewerkingen worden beiden met HIsat2 uitgevoerd.
HISAT2 is een gevoelig en snel alignment programma voor het in kaart brengen van ngs reads dit is voor zowel DNA als RNA naar een populatie van menselijke genomen en naar een enkel referentiegenoom. Gebaseerd op een uitbreiding van BWT voor grafieken (Sirén et al. 2014) is er een graph FM index (GFM) ontworpen en geïnplementeerd, een originele aanpak en de eerste implementatie ervan. Naast het gebruik van één globale GFM index die een populatie van menselijke genomen vertegenwoordigt, gebruikt HISAT2 een grote set van kleine GFM indexen die gezamelijk het gehele genoom bestrijkt. Deze kleine indexen ook wel lokale indexen genoemd, gecombindeerd met verschillende uitlijningsstragieën, maken snelle en nauwkeurige alignments van sequencing reads mogelijk. Dit nieuwe indexeringsschema wordt een Hierarchical Graph FM index (HGFM) genoemd.
#Hisat build wordt dus gebruikt om de ref genoom correct te aligneren de -p optie laat HISAT2 een specifiek aantal parallel zoek threads gebruiken die 60 zijn door ons aangegeven
hisat2-build -p 60 /students/2024-2025/Thema05/3dconformatieChromatine/Mapping_ref/ncbi_dataset/ncbi_dataset/data/GCF_000001635.27/GCF_000001635.27_GRCm39_genomic.fna hisat2_indexRead mappig:
Hisat index
Hier wordt met hisat2 de sequenties ge aligned aan het zonet gemaakte referentie genoom. De -p optie laat HISAT2 een specifiek aantal parallel zoek threads gebruiken die 80 zijn door ons aangegeven. -x geeft de basename van de index voor het reference genoom, in ons geval de hisat2_index. -1 geeft door komma’s gescheiden lijst met bestanden die mate 1’s bevatten (bestandsnaam bevat meestal _1), bijvoorbeeld -1 flyA_1.fq,flyB_1.fq. Sequenties die met deze optie worden opgegeven, moeten bestand-voor-bestand en lees-voor-lees overeenkomen met die welke zijn opgegeven in <m2>. Lezen kan een mix van verschillende lengtes zijn. Als - is opgegeven, leest hisat2 de mate 1’s van de “standaard in” of “stdin”-bestandshandle.
Daarna leest hisat de forward- en reverse paired reads de {} geeft aan dat de bestand naam alles mag zijn maar het moet _forward_paired.fastq of _rev_paired.fastq bevatten om een file te zijn die mee word genomen in de hisat.
-S zegt dat Bestand om SAM-uitlijningen naar te schrijven. Standaard worden uitlijningen naar de “standard out” of “stdout” filehandle geschreven (d.w.z. de console).
cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
hisat2 -p 80 -x hisat2_index -1 /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/trimmomatic_output/paired/{}_forward_paired.fastq -2 /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/trimmomatic_output/paired/{}_rev_paired.fastq -S /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/hisat2_mapped/{}.sam &De output van de read mapping is een sam bestand, een sam bestand ziet er zo uit:
Hier is het format te zien van een sam file. Er staan verschillende dingen in zoals de alignement info, de sequentie read, mapping quality read en nog verschillende dingen die belangrijk zijn voor de tools in de pipeline.
Sorteren en indexeren van .sam
Omzetten van .sam naar .bam: De reden om de file van een .sam naar een .bam om te zetten is omdat .bam files beter manipuleer baar zijn en ook sneller dus daarom beter zijn om mee te werken.
#Cat leest elke File-parameter in volgorde van
cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
samtools view -@40 -b /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/hisat2_mapped/{}.sam > /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/{}.bam &
Hier gebruiken we samtools view dit zet de .sam file over naar een .bam, de -@ optie word gebruikt om meerdere threads te selecteren in ons geval 40 de -b optie geeft aan welke soort file wij als output willen en -b staat natuurlijk gelijk aan .bam De andere twee lines geven de path aan waar de input files staan, en waar de output heen moet. & geeft aan dat het programma in de achtergrond moet runnen.
Index maken: Indexeer coördinaten-gesorteerde BGZIP-gecomprimeerde SAM-, BAM- of CRAM-bestanden voor snelle willekeurige toegang.Deze index is nodig wanneer regio-argumenten worden gebruikt om de samtools-weergave en soortgelijke opdrachten te beperken tot bepaalde regio’s van belang.
samtools index roept natuurlijk de index aan, de -@ optie word gebruikt om meerdere threads te selecteren in ons geval 80 en alles wordt in parallel gerund dus meerdere opdrachten tegelijk.
Verder is nog het output path te zien, hier worde
cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
samtools index -@80 /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/sorted_{}.bam
Samtools sort:
Sorteer uitlijningen op meest linkse coördinaten, op gelezen naam wanneer -n of -N worden gebruikt. De gesorteerde uitvoer wordt standaard naar de standaarduitvoer geschreven of naar het opgegeven bestand (out.bam) wanneer -o wordt gebruikt.
samtools sort roept dus de sorting aan. @ zorgt dat er een bepaald aantal threads worden gebruikt in ons geval 80. -o geeft aan dat het ouput bestand een bepaald bestand type moet zijn dit is bij ons BAM . -n geeft aan dat ons bestand op naam gelezen word.
cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
samtools sort -@80 -O BAM -n /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/{}.bam -o /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/sorted_{}.bam
samtools fixmate
samtools fixmate – vult in mate coördinaten en insert “size fields” -m voegt mate scored toe en deze zijn nodig bij markdup om de beste reads te behouden -@ 80 zegt dat er 80 threads moeten worden gebruikt voor deze opdracht verder zijn er nog 2 paths die de input en output geven.
cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
'samtools fixmate -m -@80 /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/sorted_{}.bam /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/fixed_mates_sorted_{}.bam'Sorteren op coördinaten, dit is de standaard dus hoeft er geen argument voor meegegeven worden.Vraag: Waarom zijn de gesorteerde bestanden kleiner dan de niet gesorteerde?
Sorteer uitlijningen op meest linkse coördinaten, op gelezen naam wanneer -n of -N worden gebruikt. De gesorteerde uitvoer wordt standaard naar de standaarduitvoer geschreven of naar het opgegeven bestand (out.bam) wanneer -o wordt gebruikt.
samtools sort roept dus de sorting aan. @ zorgt dat er een bepaald aantal threads worden gebruikt in ons geval 80. -o geeft aan dat het ouput bestand een bepaald bestand type moet zijn dit is bij ons BAM .
cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
'samtools sort -@80 /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/fixed_mates_sorted_{}.bam -o /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/fixed_mates_sorted_coords_{}.bam'Duplicaten verwijderen: samtools markdup markeert de duplicaten in een gecoördineerd en gesorteerde file
samtools markdup roept de markdup functie aan in samtools. -@ 80 zegt dat er 80 threads moeten worden gebruikt voor deze opdracht. -r zegt dat de gemarkeerde duplicates moeten worden verwijderd. -s zegt dat er wat basis statistieken er bij moeten worden gevoegd.
cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
'samtools markdup -@80 -r -s /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/fixed_mates_sorted_coords_{}.bam /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/dedup_{}.bam'Exploratory Data Analysis:
Featurecounts:bron: https://subread.sourceforge.net/featureCounts.html-T 64 was het maximale aantal threads dat het programma kan gebruiken.
featureCounts roept het programma feature count aan.
-T geeft de aantal threads aan die het programma mag gebruiken, in dit geval 64 dit is de maximale aantal threads dat het programma kan gebruiken door het limiet van feature counts. dit is verder te lezen in de bron hier boven.
-B creërt 1 blok index. De bouw index zal niet worden gesplitst in meerdere stukken. Hoe meer blocken een index heeft, deste langzamer de mapping snelheid. Deze optie zal de -M overschrijven als deze ook word meegegeven.
-p specifeert dat de data paired-end reads zijn, feature counts zal stoppen wanneer het een single-end reads ook wel anders dan gespecificeerd. Om de fragmenten te tellen voor paired -end reads, zal de —countReadPairs parameter moeten worden gespecificeerd
-t dit specifeert de type input van sequencing data. er zijn verschillende value’s die de volgende bevatten, 0 = RNA-seq, 1 = DNA-seq data. De gebruiker moet een bepaalde waar specificeren dit kan echter ook een karakter value zijn. Bij ons is exon aangegeven.
-g specificeert de attribute type gebruikt voor het grouperen van features (eg. exons) naar meta-features (eg. genes) wanneer GTF annotatie “provided”. gene id is de default. Dit attribute type is meestal de gen identifier. Dit argument is handig voor meta feature level samenvatting
-a geeft aan dat er een gen annotatie file is die de chorosomale coördinaten van de exonen van elk gen bevat. De path hiervan komt dan ook na de -a.
-o specifeert de base naam van de index die gemaakt word.
verder zijn dus de paths te zien.
cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
'featureCounts -T 64 -p --countReadPairs -B -t exon -g gene_id -a /students/2024-2025/Thema05/3dconformatieChromatine/snpEff/snpEff/data/mm39/genes.gtf -o /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/featurecounts_output/{}_counts.txt /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/dedup_{}.bam'
# Test:
featureCounts -T 64 -p --countReadPairs -a /students/2024-2025/Thema05/3dconformatieChromatine/snpEff/snpEff/data/mm39/genes.gtf -o /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/featurecounts_output/count.out /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/dedup_*
Resultaten van FeatureCounts:
| A s s i g n e d | 7 , 4 7 8 , 2 7 5 | 1 1 , 5 1 9 , 0 3 3 | 4 0 1 , 9 1 9 | 1 , 7 9 8 , 7 8 3 | 1 , 0 3 4 , 2 0 6 | 1 , 8 1 5 , 6 8 5 | 2 , 1 2 3 , 1 1 7 | 1 , 5 5 9 , 9 6 5 |
|---|---|---|---|---|---|---|---|---|
| U n a s s i g n e d _ U n m a p p e d |
1 0 1 , 5 6 9 | 9 1 , 9 5 7 | 3 0 3 , 1 4 7 | 5 0 1 , 9 1 8 | 1 , 9 0 0 , 7 6 5 | 2 , 0 0 6 , 3 6 6 | 2 , 3 0 5 , 5 5 2 | 3 , 1 5 3 , 7 0 0 |
| U n a s s i g n e d R e a d T y p e |
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| U n a s s i g n e d _ S i n g l e t o n |
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| U n a s s i g n e d _ M a p p i n g Q u a l i t y |
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| U n a s s i g n e d _ C h i m e r a |
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| U n a s s i g n e d _ F r a g m e n t L e n g t h |
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| U n a s s i g n e d _ D u p l i c a t e |
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| U n a s s i g n e d _ M u l t i M a p p i n g |
1 , 2 7 8 , 1 8 4 | 1 , 2 7 8 , 6 6 7 | 1 , 4 2 3 , 5 8 5 | 2 , 4 4 4 , 3 4 2 | 1 , 7 2 8 , 7 3 8 | 1 , 7 1 8 , 8 6 9 | 1 , 8 6 8 , 8 4 9 | 1 , 9 7 3 , 5 7 6 |
| U n a s s i g n e d _ S e c o n d a r y |
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| U n a s s i g n e d _ N o n S p l i t |
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| U n a s s i g n e d _ N o F e a t u r e s |
7 2 3 , 9 8 2 | 9 9 9 , 0 0 1 | 5 8 6 , 7 1 9 | 1 , 2 8 4 , 5 7 6 | 4 2 4 , 9 2 5 | 5 7 4 , 5 0 4 | 5 6 9 , 8 1 1 | 7 4 4 , 2 9 8 |
| U n a s s i g n e d O v e r l a p p i n g L e n g t h |
0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| U n a s s i g n e d _ A m b i g u i t y |
3 0 , 5 7 5 | 3 2 , 5 9 8 | 3 1 , 9 2 3 | 1 7 , 5 9 8 | 1 2 , 4 0 6 | 1 5 , 2 7 6 | 1 8 , 9 3 9 | 1 4 , 7 8 3 |
Table 2: feature count resultaten
Assigned (A): Dit zijn de reads die correct zijn toegewezen aan genen. Hoge waarden hier wijzen op een succesvolle toewijzing en een efficiënte sequencing-run.
Unassigned: Er zijn meerdere subcategorieën hier die de oorzaken van de niet-toewijzing weergeven:
Mapping Quality: Reads die door lage mapping-kwaliteit niet konden worden toegewezen.
Ambiguity: Reads die aan meerdere genen kunnen worden toegewezen en daardoor ambigu zijn.
Multi Mapping: Reads die meerdere matches hebben in het genoom.
Secondary en Duplicate: Overbodige of gedupliceerde reads, vaak door sequencing-artifacten.
No Features: Reads die niet aan een feature zijn gekoppeld, mogelijk door technische redenen.
We kunnen met de uitkomt van deze data op het eerste gezicht zien dat de Assigned waarden hoger zijn dan Unnasigned wat dus betekent dat de data redelijk goed is.
Visualisatie van featurecounts
Het van FeatureCounts verkregen .csv bestand kan nu als tabel in R geladen worden. De namen van de samples zijn erg lang en bevatten het gehele pad naar de bestanden, daarom worden de kolomen van de samples met de onderstaande code hernoemd naar de SSR namen.
#counts <- read.delim(file = "/home/floris/Documenten/Data_set/PRJNA885415/Transcriptomics/featurecounts/featurecounts_output.csv", header = F)
counts = read.csv( "/Users/jarnoduiker/Desktop/Transcriptomics_data/featurecounts_output.csv", row.names = "Geneid", sep="", head=T, skip=1,)
# Geef namen aan de kolomen
names(counts) <- c("Chr","Start","End","Strand","Length","SRR21754408","SRR21754417","SRR21754418","SRR21754419","SRR21754420","SRR21754421","SRR21754422","SRR21754423")
# conflict_prefer("select", "dplyr") is nodig omdat er blijkbaar twee libraries zijn die de functie "select" aanbieden.
conflict_prefer("select", "dplyr")## [conflicted] Will prefer dplyr::select over any other package.counts_df <- counts %>%
select(-c(Chr, Start, End, Strand, Length))
glimpse(counts_df)## Rows: 41,079
## Columns: 8
## $ SRR21754408 <int> 0, 0, 1, 0, 10, 1, 0, 0, 0, 1, 438, 862, 1, 839, 10, 0, 43…
## $ SRR21754417 <int> 0, 0, 1, 1, 10, 0, 0, 0, 0, 0, 606, 1217, 3, 1244, 23, 0, …
## $ SRR21754418 <int> 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 14, 20, 0, 36, 1, 0, 9, 0, 0…
## $ SRR21754419 <int> 0, 0, 1, 0, 2, 0, 0, 0, 0, 0, 66, 127, 1, 140, 3, 0, 80, 0…
## $ SRR21754420 <int> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 44, 114, 0, 117, 0, 0, 41, 0…
## $ SRR21754421 <int> 0, 0, 1, 3, 1, 0, 1, 0, 1, 0, 61, 151, 0, 213, 5, 1, 85, 0…
## $ SRR21754422 <int> 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 96, 209, 0, 240, 1, 0, 131, …
## $ SRR21754423 <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 60, 144, 0, 141, 1, 0, 94, 0…Hier wordt de colname veranderd naar groep, en de coldata word een data frame met de aangegeven samples, (SRR21754408,SRR21754417,SRR21754418,SRR21754419,SRR21754420,SRR21754421,SRR21754422,SRR21754423) daarna worden de colnames weer aan gepast naar groep in het data frame en de kollom groep wordt ook nog factor gebruikt om een “vector” naar factor te transformen.
colnames <- c("Groep")
coldata <- data.frame(c("SRR21754408","SRR21754417","SRR21754418","SRR21754419","SRR21754420","SRR21754421","SRR21754422","SRR21754423"), c("Jong","Jong","Jong","Jong","Oud","Oud","Oud","Oud"), row.names = 1)
colnames(coldata) <- colnames
#coldata <- coldata[,c("Groep")]
coldata$Groep <- factor(coldata$Groep)
#coldata$Type <- factor(coldata$Type)Hier word er door het data frame geloopt en met pivot longer dit “pivot” een data frame van wijdformaat naar langformaat.
Dan wordt er ggplot gebruikt om een box plot te maken om te laten zien hoeveel reads er bij elke sample zit.
counts_long <- counts_df %>%
pivot_longer(cols = everything(),
names_to = "sample_naam",
values_to = "reads")
ggplot(data = counts_long, aes(x = sample_naam, y = reads)) +
geom_boxplot() +
xlab(" ") +
ylab("Aantal reads:") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
Figuur 8: boxplot reads per sample
In het plot is te zien dat het grootste deel van de samples tussen de 0 en 50000 reads zit echter is er bij samples 4419 te zien dat er 1 uitschieter is van boven de 2000000 wat wel interessant is want waarom zitten hier zoveel reads? Het kan meerdere dingen zijn van technische fouten tot biologische variant.
Echter ziet het er verder redelijk goed uit waardoor we verder kunenen
In het volgende codeblok worden de kolommen geselecteerd waarin de jonge en oude muizen staan. Deze worden in nieuwe variabellen gezet zodat we ze afhankelijk van elkaar kunnen bekijken en dus makkelijker vergelijken.
oude_muizen <- counts_df[,c(5:8)]
jonge_muizen <- counts_df[,c(1:4)]Hier word de data van elke sample die hoort bij jong en oud in hun toebehorende data frame gezet.
# DF met alle rijen waarbij de jonge muizen een coverage hoger dan 0 hebben voor alle samples:
jonge_muizen_df <- counts[which(counts$SRR21754408 > 0 & counts$SRR21754418 > 0 & counts$SRR21754419 & counts$SRR21754417), ]
# DF met alle rijen waarbij de oude muizen een coverage hoger dan 0 hebben voor alle samples:
oude_muizen_df <- counts[which(counts$SRR21754420 > 0 & counts$SRR21754421 > 0 & counts$SRR21754422 & counts$SRR21754423 > 0), ]In het volgende stuk code is een grafiek weergegeven van alle genen met een hogere expressie dan 0.
#ggplot word hier aangeroepen met ggplot(), dan is jonge_muizen_df als de data source gebruikt. in de aes moeten de x en y worden gespecifeerd waar de x row.names is wat dus de gen namen zijn en y de lengte is wat dus de groote is van het gen. daarna is geom_bar() aangeroepen dit geeft aan dat de plot een barplot moet zijn. xlab en ylab verzorgen de tekst op de x en y as. daarna nog het thema van hoe de grafiek er moet uitzien
ggplot(data = jonge_muizen_df, mapping = aes(x = row.names(jonge_muizen_df) , y = Length)) +
geom_bar(stat = "identity") +
xlab(" ") +
ylab("Gene length: (bp)") +
theme(axis.text.x = element_blank())
Figuur 9: gen expressie jonge muizen
Deze plot laat de lengte van genen in baseparen (bp) zien (zoals te zien op de y as), elke balk vertegenwoordigt een gen en geeft de lengte aan. op de x-as staan de individuele genen alleen als deze weergegeven worden word de grafiek onleesbaar. Wat opvalt is dat de meeste genen relatief kort zijn vaak onder de 50.000, maar er zijn een paar genen die veel langer zijn en uitschieters vormen deze kunnen belangrijke / complexe genen zijn die vaak meer exons en introns bevatten of belangrijke regulatorische sequenties dit is dus alleen in jonge muizen.
#ggplot word hier aangeroepen met ggplot(), dan is jonge_muizen_df als de data source gebruikt. in de aes moeten de x en y worden gespecifeerd waar de x row.names is wat dus de gen namen zijn en y de lengte is wat dus de groote is van het gen. daarna is geom_bar() aangeroepen dit geeft aan dat de plot een barplot moet zijn. xlab en ylab verzorgen de tekst op de x en y as. daarna nog het thema van hoe de grafiek er moet uitzien
ggplot(data = oude_muizen_df, mapping = aes(x = row.names(oude_muizen_df) , y = Length)) +
geom_bar(stat = "identity") +
xlab(" ") +
ylab("Gene length: (bp)") +
theme(axis.text.x = element_blank())
Figuur 10: gen expressie oude muizen
Deze plot laat de lengte van genen in baseparen (bp) zien (zoals te zien op de y as), elke balk vertegenwoordigt een gen en geeft de lengte aan. op de x-as staan de individuele genen alleen als deze weergegeven worden word de grafiek onleesbaar. net zoals de vorige plot echter is wel te zien dat de lengte van de genen hier wel anders zijn dan die bij de jonge muizen dit is te zien aan de rechter kant van de x-as hier zijn 2 grote pieken te zien die niet in de jonge grafiek is te zien. Dit kan wijzen op verandering in gen lengte - sequentie lengte met leeftijd.
myColors <- hue_pal()(4)
## Plot the log2-transformed data with a 0.1 pseudocount
plotDensity(log2(counts_df + 0.1), col=rep(myColors, each=3),
lty=c(1:ncol(counts_df)), xlab='Log2(count)',
main='Expression Distribution')
## Add a legend and vertical line
legend('topright', names(counts_df), lty=c(1:ncol(counts_df)),
col=rep(myColors, each=3))
abline(v=-1.5, lwd=1, col='red', lty=2)
Figuur 11: expressie distributie
Deze plot laat de expressie verdeling van verschillende samples (SRR21754408, SRR21754417, enzovoort) zien met log2count. Wat de log-getransformeerde read counts per gen weergeeft. op de x-as staat de genexpressie weergegeven als log2count, terwijl de y-as de dichtheid aangeeft.
de pieken rond de negatieve waarden geven waarschijlijk genen aan met lage of geen expressie in de meeste samples waar de Log2(count)-waarden positief zijn, worden genen met hogere expressie zichtbaar, hoewel dit minder frequent voorkomt.
Wat is uit te laten is dat er een grote negatieve piek is. Een negatieve log2 zegt dat er negatieve of geen expressie. De meeste genen die in onze data zitten hebben dus minder of geen expressie en zijn dus niet super interessant.
Verder is een kleine piek te zien bij 0 dit zijn dus genen die in DEZE samples niet expressed worden. als laatste is te zien dat de 4408 en 4417 sample bij de 10 ook nog expressie hebben dus bij deze samples zouden er genen kunnen zijn die overexpressed zijn.
Het aantal reads met een coverage van meer dan 0 van de jonge muizen: 14315
length(which(counts$SRR21754408 > 0 & counts$SRR21754418 > 0 & counts$SRR21754419 & counts$SRR21754417 > 0))## [1] 14290Het aantal reads met een coverage van meer dan 0 van de oude muizen: 15247
length(which(counts$SRR21754420 > 0 & counts$SRR21754421 > 0 & counts$SRR21754422 & counts$SRR21754423 > 0))## [1] 15247deseq2 Inlezen van de feature counts data inclusief de metadata collomen De deseq is gemaakt door Ivar in zijn logboek zal er uitgebreide functie uitleg zitten. DE DESEQ DATA IS NIET HELEMAAL CORRECT DIT IS OMDAT ER 1 SAMPLE MIST: SRR21754417
Deze deseq2 is gemaakt door ivar. om hier betere en verdere uitleg van te krijgen kijk naar het logboek van ivar.
# inlezen
# skip 2 lines met import data
transcrip_table <- read.delim("/Users/jarnoduiker/Desktop/Transcriptomics_data/count.out",header = FALSE,skip = 2)
# names vector van floris
names(transcrip_table) <- c("Gene_id","Chr","Start","End","Strand","Length","SRR21754408","SRR21754418","SRR21754419","SRR21754420","SRR21754421","SRR21754422","SRR21754423")
#print
head(transcrip_table)## Gene_id
## 1 Gm26206
## 2 Xkr4
## 3 Gm53491
## 4 Rp1
## 5 Sox17
## 6 Gm22307
## Chr
## 1 NC_000067.7
## 2 NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7
## 3 NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7
## 4 NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7
## 5 NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7
## 6 NC_000067.7
## Start
## 1 3172239
## 2 3269956;3283662;3284705;3393039;3420932;3491925;3491925;3491925;3491925;3740775;3740775;3740775;3740775
## 3 3363896;3365064;3396936;3422903;3436480;3438011;3438565;3438565
## 4 4185896;4189891;4190089;4190238;4212835;4212835;4218035;4218035;4218835;4218835;4234078;4234078;4240428;4240428;4267757;4267757;4276883;4276883;4296834;4296834;4298666;4298666;4301276;4301276;4313356;4313356;4313640;4313640;4313766;4313766;4315254;4315254;4331750;4331750;4337692;4337692;4354989;4354989;4361069;4363149;4363149;4381493;4381493;4414373;4422133;4422133;4422133;4422133;4422425;4422425;4422425;4422425;4430423;4430423;4479393;4479393
## 5 4561154;4561154;4561154;4561154;4561154;4561154;4561154;4561154;4561154;4561154;4563323;4563323;4563323;4563323;4563323;4563323;4563995;4563995;4563995;4563995;4563995;4563995;4565359;4565359;4565359;4565359;4565359;4566514;4566514;4566514;4566514;4566514;4566514;4566514;4566514
## 6 4599240
## End
## 1 3172348
## 2 3277540;3287191;3287191;3393983;3422536;3492124;3492124;3492124;3492124;3741733;3741733;3741733;3741721
## 3 3364050;3365091;3397056;3438060;3436728;3438060;3448035;3448035
## 4 4186171;4189935;4190296;4190296;4212989;4212989;4218186;4218186;4218967;4218967;4234164;4234164;4240627;4240627;4267864;4267864;4277060;4277060;4297046;4297046;4298842;4298842;4301367;4301367;4313485;4313485;4313671;4313671;4313842;4313842;4315329;4315329;4331828;4331828;4337843;4337843;4355121;4355121;4363235;4363235;4363235;4381656;4381656;4420314;4422304;4422304;4422304;4422304;4423060;4423060;4423060;4423060;4430526;4430526;4479508;4479489
## 5 4562891;4562891;4562891;4562891;4562891;4562891;4562891;4562891;4562891;4562891;4563831;4565421;4563689;4563689;4563713;4565617;4564086;4564086;4564086;4564086;4564086;4564086;4565617;4566165;4566165;4565421;4565421;4566596;4566664;4566664;4566664;4566664;4566664;4567624;4566664
## 6 4599346
## Strand
## 1 +
## 2 -;-;-;-;-;-;-;-;-;-;-;-;-
## 3 +;+;+;+;+;+;+;+
## 4 -;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-
## 5 -;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-
## 6 +
## Length SRR21754408 SRR21754418 SRR21754419 SRR21754420 SRR21754421
## 1 110 0 0 0 0 0
## 2 14824 0 0 1 0 0
## 3 24933 1 1 1 1 1
## 4 12004 0 1 0 0 0
## 5 5692 10 10 0 2 0
## 6 107 1 0 0 0 0
## SRR21754422 SRR21754423 NA
## 1 0 0 0
## 2 0 1 0
## 3 1 1 0
## 4 3 0 0
## 5 1 1 0
## 6 0 1 0De import data veranderen naar een count matrix voor de deseq functie
Er zijn hier een paar dingen die gebeuren. Er is een nieuwe variable count_matrix, hierin word een data frame gezet die ook gelijk word aangemaakt. in dit data frame zit transcrip_table id en alle samples. De names worden ook veranderd zodat het goed leesbaar is.
count_matrix <- data.frame(transcrip_table$Gene_id, transcrip_table$SRR21754408,transcrip_table$SRR21754418,transcrip_table$SRR21754419,transcrip_table$SRR21754420,transcrip_table$SRR21754421,transcrip_table$SRR21754422,transcrip_table$SRR21754423)
names(count_matrix) <- c("Gene_ID","SRR21754408","SRR21754418","SRR21754419","SRR21754420","SRR21754421","SRR21754422","SRR21754423")
head(count_matrix)## Gene_ID SRR21754408 SRR21754418 SRR21754419 SRR21754420 SRR21754421
## 1 Gm26206 0 0 0 0 0
## 2 Xkr4 0 0 1 0 0
## 3 Gm53491 1 1 1 1 1
## 4 Rp1 0 1 0 0 0
## 5 Sox17 10 10 0 2 0
## 6 Gm22307 1 0 0 0 0
## SRR21754422 SRR21754423
## 1 0 0
## 2 0 1
## 3 1 1
## 4 3 0
## 5 1 1
## 6 0 1Hier wordt weer een dataframe aangemaakt voor de nieuwe variabele conditie_matrix. hierin staat de sample, jong of oud en de PE. Daarna word het data frame getransformeerd door de sample names niet als kol naam te zetten maar als row names en de col names aan passen naar sample, group en type. verder word de data omgezet naar “factor” van int / chr
# condition data frame
conditie_matrix <- data.frame(
c("SRR21754408","young","PE"),
c("SRR21754418","young","PE"),
c("SRR21754419","young","PE"),
c("SRR21754420","old","PE"),
c("SRR21754421","old","PE"),
c("SRR21754422","old","PE"),
c("SRR21754423","old","PE")
)
transpose <- data.frame(t(conditie_matrix))
# corigeer frame voor input functie deseq2
# condities moeten factoren zijn voor input functie
conditie_matrix <- data.frame(as.factor(transpose$X1),as.factor(transpose$X2),as.factor(transpose$X3))
names(conditie_matrix) <- c("sample","group","type")
head(conditie_matrix)## sample group type
## 1 SRR21754408 young PE
## 2 SRR21754418 young PE
## 3 SRR21754419 young PE
## 4 SRR21754420 old PE
## 5 SRR21754421 old PE
## 6 SRR21754422 old PEHier wordt het cts data frame aangemaakt waar alle data van de samples ingezet word, hierin staat eigenlijk of de gen aanwezig is in de sample
# samples in een nieuw dataframe
cts <- data.frame(count_matrix$SRR21754408, count_matrix$SRR21754418, count_matrix$SRR21754419, count_matrix$SRR21754420, count_matrix$SRR21754421, count_matrix$SRR21754422, count_matrix$SRR21754423)
# verander colnames en rownames
# colnames moet gelijk zijn aan rownames conditie matrix
colnames(cts) <- c("SRR21754408","SRR21754418","SRR21754419","SRR21754420","SRR21754421","SRR21754422","SRR21754423")
# rownames voor de gene names
row.names(cts) <- count_matrix$Gene_ID
head(cts)## SRR21754408 SRR21754418 SRR21754419 SRR21754420 SRR21754421 SRR21754422
## Gm26206 0 0 0 0 0 0
## Xkr4 0 0 1 0 0 0
## Gm53491 1 1 1 1 1 1
## Rp1 0 1 0 0 0 3
## Sox17 10 10 0 2 0 1
## Gm22307 1 0 0 0 0 0
## SRR21754423
## Gm26206 0
## Xkr4 1
## Gm53491 1
## Rp1 0
## Sox17 1
## Gm22307 1Al deze voorgemaakte data frames komen nu samen om het deseq object te maken.
Deseq object aanmaak
# Cts de verwerkte fetuere counts matrix
# conditie matrix is de metadata matrix voor de vergelijking
# design is de colom (datatype=factor) van de conditie matrix.
dds <- DESeqDataSetFromMatrix(countData = cts,
colData = conditie_matrix,
design = ~ group)
dds## class: DESeqDataSet
## dim: 41079 7
## metadata(1): version
## assays(1): counts
## rownames(41079): Gm26206 Xkr4 ... TrnT TrnP
## rowData names(0):
## colnames(7): SRR21754408 SRR21754418 ... SRR21754422 SRR21754423
## colData names(3): sample group typedeseq analyse object
deseq <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingresults(deseq)## log2 fold change (MLE): group young vs old
## Wald test p-value: group young vs old
## DataFrame with 41079 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## Gm26206 0.000000 NA NA NA NA NA
## Xkr4 0.986635 0.460688 2.72789 0.168881 0.865890 NA
## Gm53491 1.726109 -1.150346 1.67134 -0.688278 0.491278 NA
## Rp1 0.533061 -1.969095 3.02457 -0.651033 0.515025 NA
## Sox17 1.314470 0.533431 1.62886 0.327486 0.743300 NA
## ... ... ... ... ... ... ...
## ND6 1373.89862 -0.470926 0.536920 -0.877088 0.3804389 0.923757
## TrnE 3.49194 0.338739 1.133636 0.298808 0.7650869 NA
## CYTB 4964.47156 -0.586623 0.279874 -2.096027 0.0360798 0.415000
## TrnT 2.75378 -0.902012 1.144438 -0.788170 0.4305971 NA
## TrnP 91.48406 -0.436424 0.337111 -1.294602 0.1954574 0.795742De resultaten van de deseq geeft collommen met verschillende info
baseMean - het gemiddelde van de genormaliseerde telwaarden, gedeeld door groottefactoren, genomen over alle monsters
Log2FoldChange - De verhouding van de genormaliseerde gemiddelde gen-UMI-tellingen in elke cluster/groep ten opzichte van alle andere clusters/groepen ter vergelijking
lfcSE - De standaardfoutschatting voor de log2-voudige veranderingsschatting
stat - De waarde van de test statistiek voor het gen of het transcript
pvalue - is dus de pvalue
padj - is de adjusted pvalue voor meerdere testen voor genen of transcripten
Wat kan je dan met deze dingen. Je kan hier verschillende dingen mee zoals een kegg pathway analyse, volcano plots en andere visualisaties. Het wordt dus gebruikt om te laten zien hoe de genen worden expressed en daarmee is te zien of er een verkeerd werkend gen is.
result deseq object waar gegroupeert word waar dus de resultaten van de oude muizen samples worden vergeleken met de jonge muizen samples.
# vergelijk groep
resultaat_deseq <- results(deseq, contrast = c("group", "young","old"))
resultaat_deseq## log2 fold change (MLE): group young vs old
## Wald test p-value: group young vs old
## DataFrame with 41079 rows and 6 columns
## baseMean log2FoldChange lfcSE stat pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
## Gm26206 0.000000 NA NA NA NA NA
## Xkr4 0.986635 0.460688 2.72789 0.168881 0.865890 NA
## Gm53491 1.726109 -1.150346 1.67134 -0.688278 0.491278 NA
## Rp1 0.533061 -1.969095 3.02457 -0.651033 0.515025 NA
## Sox17 1.314470 0.533431 1.62886 0.327486 0.743300 NA
## ... ... ... ... ... ... ...
## ND6 1373.89862 -0.470926 0.536920 -0.877088 0.3804389 0.923757
## TrnE 3.49194 0.338739 1.133636 0.298808 0.7650869 NA
## CYTB 4964.47156 -0.586623 0.279874 -2.096027 0.0360798 0.415000
## TrnT 2.75378 -0.902012 1.144438 -0.788170 0.4305971 NA
## TrnP 91.48406 -0.436424 0.337111 -1.294602 0.1954574 0.795742Het is moeilijk uit te lezen op deze manier dus we zetten het in een variable zodat het kan vergeleken worden met latere visualisaties.
LFC shrinkage object
lfcShrink - gebruikt om de log fold changes (LFC) van genexpressieniveaus te verkleinen of “in te schatten”. Dit komt vaak voor in analyses met DESeq2. dit doe je omdat na de standaard DEseq2-analyse sommige genen extreem hoge of lage LFC-waarden kunnen hebben. vooral als deze lage reads hebben of grote variaties. Zulke extreme waarden kunnen ruis toevoegen aan de resultaten en om dit te vermijden past lcfShrink shrinkage toe op de LFC waarden wat dus betekent dat het deze waarden dichter naar 0 toe trekt voor genen met lage betrouwbaarheid (de genen met weinig reads of hoge variantie) dit maakt de resultaten stabieler en helpt focus op de genen met een goed-geschatte LFC-waarden
deseq_shrink <- lfcShrink(deseq, coef = "group_young_vs_old", type = "apeglm")## using 'apeglm' for LFC shrinkage. If used in published research, please cite:
## Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
## sequence count data: removing the noise and preserving large differences.
## Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895deseq_shrink## log2 fold change (MAP): group young vs old
## Wald test p-value: group young vs old
## DataFrame with 41079 rows and 5 columns
## baseMean log2FoldChange lfcSE pvalue padj
## <numeric> <numeric> <numeric> <numeric> <numeric>
## Gm26206 0.000000 NA NA NA NA
## Xkr4 0.986635 0.00282996 0.127617 0.865890 NA
## Gm53491 1.726109 -0.00384646 0.127515 0.491278 NA
## Rp1 0.533061 -0.00356699 0.127684 0.515025 NA
## Sox17 1.314470 0.00306657 0.127293 0.743300 NA
## ... ... ... ... ... ...
## ND6 1373.89862 -0.02464090 0.127761 0.3804389 0.923757
## TrnE 3.49194 0.00457367 0.126905 0.7650869 NA
## CYTB 4964.47156 -0.14708214 0.227104 0.0360798 0.415000
## TrnT 2.75378 -0.01102058 0.127642 0.4305971 NA
## TrnP 91.48406 -0.05977954 0.137425 0.1954574 0.795742hoeveelheid p waardes onder de 0.01 van de 41079
hier uit is te halen dat de 43 getelde sequenties zeer hoge statistische significantie hebben.
sum(resultaat_deseq$padj < 0.01, na.rm=TRUE)## [1] 43Resultaat p waardes onder de 0.05
Hier worden de resultaten van deseq gefilterd op een p-waarde van < 0.05 zodat deze later kunnen worden gebruikt.
deseq_05 <- results(deseq, alpha=0.05)
summary(deseq_05)##
## out of 28289 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 47, 0.17%
## LFC < 0 (down) : 138, 0.49%
## outliers [1] : 7, 0.025%
## low counts [2] : 16062, 57%
## (mean count < 9)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?resultsma plot pre shrinkage
plotMA(resultaat_deseq, ylim=c(-2,2))
Figuur 12: MA plot eerste deseq zonder shrinkage
ma plot met shrinkage
plotMA(deseq_shrink, ylim=c(-3,3))
Figuur 13: Ma plot eerste deseq met shrinkage
Een MA-plot is een plot die de verhouding tussen genexpressie in twee condities laat zien. Het helpt om te identificeren welke genen significant differentieel tot expressie komen. De genen die gekleurd zijn interessant, dit is omdat ze significante genexpressie hebben dus meer expressie dan de norm (de grijze bolletjes en lijn)
volcano plot
Een volcano plot is een visualisatietool die genen laat zien die significant en sterk differentieel tot expressie komen tussen twee condities. De grootte van het expressieverschil (fold change) en de significantieniveau (padj) van elk gen, waardoor je snel kunt zien welke genen biologisch relevant zijn.
EnhancedVolcano(deseq_shrink, x="log2FoldChange", y = "padj",lab = rownames(deseq_shrink), title = "young versus old", pCutoff = 0.05,FCcutoff = 2, ylim = c(0,50),legendPosition = "none",subtitle = "FDR <= 0.05 and absolute FC >= 2")
Figuur 14:Volcano plot jong vs oud
Hier is onze volcano plot, er zijn een paar dingen te zien, het grootste deel van de genen zit tussen de standaard opgestelde voorwaarden voor singificantie (de 3 onderbroken lijnen) degene die er niet tussen zitten zijn rood gekleurd. Er zijn een paar namen te zien zoals Klra8 en Clnk wat niet de pax5 en ebf1 zijn die wij onderzoeken.
plotCounts(dds, "Pax5", intgroup = "group")
Figuur 15: pax5 expressie
Hier is de genexpressie van de pax5 gen te zien per groep. De expressie van de pax5 gen is veel hoger gemiddeld in de jongere muizen wat verwacht is omdat jongere muizen nog veel ontwikkelen. er is wel 1 uitschieter bij de oudere muizen van 1600+ counts. Dit is interessant echter zou ik dit niet nu kunnen verklaren omdat er weinig verdere info is over de muizen behalve waar ze vandaan komen en dus kan het een uitschieter zijn zonder echte waarde.
plotCounts(dds, "Ebf1", intgroup = "group")
Figuur 16: expressie Ebf1 jong en oud
Hier is de gen epxressie te zien van de Ebf1 De expressie van de Ebf1 gen is veel hoger gemiddeld in de jongere muizen wat verwacht is omdat jongere muizen nog veel ontwikkelen. er is weer 1 uitschieter bij de oudere muizen van 6000-7000 counts wat waarschijnlijk dezelfde muis is als de vorige uitschieter en dit is dus niet uit te leggen op dit moment.
exctracting transformd value voor pca
vsd <- vst(dds, blind = FALSE)
head(assay(vsd))## SRR21754408 SRR21754418 SRR21754419 SRR21754420 SRR21754421 SRR21754422
## Gm26206 5.033995 5.033995 5.033995 5.033995 5.033995 5.033995
## Xkr4 5.033995 5.033995 5.641393 5.033995 5.033995 5.033995
## Gm53491 5.161996 5.136355 5.641393 5.330047 5.405110 5.308396
## Rp1 5.033995 5.136355 5.033995 5.033995 5.033995 5.507854
## Sox17 5.437583 5.357079 5.033995 5.451948 5.033995 5.308396
## Gm22307 5.161996 5.033995 5.033995 5.033995 5.033995 5.033995
## SRR21754423
## Gm26206 5.033995
## Xkr4 5.287323
## Gm53491 5.287323
## Rp1 5.033995
## Sox17 5.287323
## Gm22307 5.287323PCA plot waarbij je kunt zien dat de groepen verschillen en dat bij jong 1 outlyer is ssr 21754408 (gevonden met het gecomente pca plot hieronder)
# standaard pca
plotPCA(vsd, intgroup="group", ntop = 40000)## using ntop=40000 top features by variance
Figuur 17: verschillen in groepen
# sample identificatie
#plotPCA(vsd, intgroup=c("group","sample"))22/10/2024 KEGG pathway
Jarno’s versie van de KEGG pathway - problemen eind visualisatie.
library(clusterProfiler)## ## clusterProfiler v4.12.6 Learn more at https://yulab-smu.top/contribution-knowledge-mining/
##
## Please cite:
##
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan,
## X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal
## enrichment tool for interpreting omics data. The Innovation. 2021,
## 2(3):100141library(org.Mm.eg.db)## Loading required package: AnnotationDbi##
## Attaching package: 'AnnotationDbi'## The following object is masked from 'package:clusterProfiler':
##
## select## The following object is masked from 'package:dplyr':
##
## select## library(DOSE)## DOSE v3.30.5 For help: https://yulab-smu.top/biomedical-knowledge-mining-book/
##
## If you use DOSE in published research, please cite:
## Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an R/Bioconductor package for Disease Ontology Semantic and Enrichment analysis. Bioinformatics 2015, 31(4):608-609#een lijst met alle gen entrez id's die wij nodig hebben is de 17681
search_kegg_organism('mmu', by='kegg_code')## kegg_code scientific_name common_name
## 26 mmur Microcebus murinus gray mouse lemur
## 30 mmu Mus musculus house mouse
## 6386 mmuc Mycolicibacterium mucogenicum Mycolicibacterium mucogenicum# hier is een eerste poging op het filteren van de data op de juiste pvalue deze was echter later veranderd omdat het gehele dataframe nodig was van de deseq resultaten
padj <- resultaat_deseq$padj
padj[is.na(padj)] <- 1
filtered_deseq_result <- rownames(resultaat_deseq)[padj < 0.05]
#head(filtered_deseq_result)Hier worden de resultaten uit de deseq gefilterd. De manier dat we filteren is door een subset te maken. deze subset moet een P adjusted value hebben van minder dan 0.05. De reden dat we voor 0.05 hebben gekozen is omdat kort gezegt Fisher suggested is en ook dus de industry standaard.
subset_resultaat_deseq <- subset(resultaat_deseq, padj <= 0.05)hier wordt de subset omgezet naar een dataframe dit is zodat er beter mee is te werken.
resultaat_deseq_datafr <- as.data.frame(subset_resultaat_deseq)
resultaat_deseq_datafr## baseMean log2FoldChange lfcSE stat pvalue
## Gm4956 9.943480 -3.3573315 0.8726760 -3.847168 1.194909e-04
## Gm41914 34.779063 -1.4053941 0.3906152 -3.597900 3.207974e-04
## Actr3 1337.546469 -0.6571379 0.1918505 -3.425261 6.142093e-04
## Rab7b 24.367399 -1.7213192 0.4883002 -3.525125 4.232827e-04
## Fasl 8.598747 -3.3404567 0.8029135 -4.160419 3.176642e-05
## Xcl1 21.771844 -2.6391352 0.6993227 -3.773845 1.607510e-04
## Fcgr2b 190.309764 -0.9868967 0.2799897 -3.524761 4.238651e-04
## Fcgr4 72.033447 -1.2747637 0.3755859 -3.394067 6.886284e-04
## Pcp4l1 94.026296 0.9150324 0.2572807 3.556553 3.757522e-04
## Slamf7 559.786247 0.7016975 0.2054318 3.415720 6.361365e-04
## Atf3 82.554561 1.6651756 0.4186536 3.977455 6.965692e-05
## Olfm1 27.374401 -1.7838464 0.4300855 -4.147656 3.358971e-05
## Gm39861 17.574010 2.9629182 0.6406062 4.625179 3.742750e-06
## Rasgrp1 30.901724 -2.0215204 0.5151597 -3.924066 8.706704e-05
## Gatm 91.795531 -1.5794845 0.4351162 -3.630029 2.833893e-04
## A630026N12Rik 242.928780 0.9573570 0.2432206 3.936168 8.279316e-05
## Dusp2 240.571632 1.5378327 0.3853270 3.990981 6.580046e-05
## Tgif2 33.922861 1.5549204 0.3881685 4.005787 6.181126e-05
## Samhd1 332.449227 -0.7677018 0.1847198 -4.156034 3.238191e-05
## Pltp 43.673993 -1.5016411 0.4239877 -3.541709 3.975437e-04
## Zfp971 114.371102 -0.9442563 0.2387260 -3.955398 7.640722e-05
## Sirpb1a 23.064096 -1.6362921 0.4606710 -3.551976 3.823502e-04
## Zbtb7b 130.096675 -1.1371122 0.2840976 -4.002541 6.266589e-05
## S100a8 5291.848808 -1.3551086 0.3066503 -4.419068 9.912751e-06
## Chil3 2856.832134 -0.9274261 0.2613861 -3.548108 3.880089e-04
## S1pr1 68.469060 0.9960447 0.2822329 3.529159 4.168829e-04
## Cnn3 242.082189 0.7292388 0.2107532 3.460155 5.398641e-04
## Gem 45.303741 1.7063905 0.4627709 3.687333 2.266165e-04
## Atp6v0d2 62.772306 -4.0316520 1.1061661 -3.644708 2.676958e-04
## Klf4 117.552904 1.2566871 0.3213334 3.910851 9.197158e-05
## Slc31a2 112.966731 -0.8872508 0.2513080 -3.530532 4.147250e-04
## Megf9 219.205452 -0.7698789 0.2237705 -3.440484 5.806743e-04
## Jun 1227.719677 0.9794741 0.2296260 4.265520 1.994374e-05
## Pde4b 470.847616 0.8981810 0.2020846 4.444579 8.806402e-06
## Zc3h12a 112.832177 0.8806101 0.2382704 3.695843 2.191585e-04
## Man1c1 72.395652 -1.2627813 0.3654650 -3.455273 5.497363e-04
## Arhgef10l 26.712355 -1.6605640 0.4875159 -3.406174 6.588010e-04
## Clnk 15.059805 -2.7826646 0.6056346 -4.594626 4.335253e-06
## Prom1 91.475427 -1.0718175 0.3158507 -3.393431 6.902303e-04
## Anxa3 163.928834 -1.0527586 0.2881733 -3.653214 2.589781e-04
## Hpse 49.523869 -1.9132288 0.4340713 -4.407638 1.045042e-05
## Btbd8 77.612103 1.2308148 0.3241262 3.797332 1.462616e-04
## Oasl2 101.904126 -0.9351533 0.2738361 -3.415011 6.377934e-04
## Ncf1 944.231018 -0.8962044 0.2517757 -3.559535 3.715115e-04
## Mgam 143.380236 -0.9896614 0.2448906 -4.041238 5.316971e-05
## Gimap3 56.929469 -1.3256432 0.3483551 -3.805436 1.415544e-04
## Gpnmb 374.813201 -6.5355986 1.2766037 -5.119520 3.063139e-07
## Skap2 153.783791 -0.7211823 0.2125233 -3.393427 6.902393e-04
## Gm44008 29.837504 1.5470014 0.4465587 3.464274 5.316653e-04
## Anxa4 62.430947 -1.5159452 0.4448542 -3.407735 6.550451e-04
## Rassf4 192.962404 -1.0461482 0.2759422 -3.791186 1.499294e-04
## C3ar1 43.136920 -1.5170326 0.4384632 -3.459886 5.404051e-04
## Clec4a3 74.859998 -1.5383463 0.4127501 -3.727065 1.937225e-04
## Klrb1a 10.476746 -2.4971028 0.7159015 -3.488054 4.865505e-04
## Cd69 1113.109624 1.4464319 0.2958760 4.888643 1.015336e-06
## Clec7a 257.033375 -0.9138534 0.2692232 -3.394408 6.877725e-04
## Klre1 15.172128 -2.2812471 0.6350108 -3.592454 3.275784e-04
## Klrk1 57.751450 -2.3554778 0.4835277 -4.871443 1.107862e-06
## Klrc3 8.230109 -3.2380495 0.9087374 -3.563240 3.663056e-04
## Klri2 29.069750 -2.6545131 0.5874070 -4.519036 6.212194e-06
## Klra8 19.849416 -3.6380844 0.7124949 -5.106120 3.288405e-07
## Klra7 19.080002 -2.8744625 0.6441256 -4.462581 8.097846e-06
## Pirb 315.668459 -1.0593741 0.2349893 -4.508181 6.538590e-06
## Ncr1 55.410115 -2.6790885 0.5715374 -4.687512 2.765462e-06
## Pglyrp1 604.584957 -1.2253921 0.3135151 -3.908559 9.284828e-05
## Cd177 1560.268810 -1.0378258 0.2933935 -3.537317 4.042143e-04
## Zfp36 901.578436 1.0917871 0.2217964 4.922475 8.545669e-07
## Rinl 125.749178 -0.9031073 0.2419821 -3.732125 1.898710e-04
## Ppp1r15a 238.156393 1.4704520 0.3159321 4.654330 3.250359e-06
## Svip 41.727998 -1.3100791 0.3823318 -3.426550 6.113001e-04
## Mfge8 92.586206 -1.7886042 0.4393962 -4.070596 4.689302e-05
## Anpep 58.089956 -3.4710640 0.9465854 -3.666932 2.454782e-04
## Gm15501 34.107488 -1.9321662 0.4947411 -3.905409 9.406617e-05
## Tenm4 36.558942 -1.8938937 0.5099552 -3.713843 2.041352e-04
## Bmal1 422.163771 0.7897508 0.2312758 3.414757 6.383891e-04
## Itpripl2 172.168920 -1.0845697 0.2924845 -3.708128 2.087975e-04
## Igsf6 232.099997 -1.1001642 0.2868806 -3.834921 1.256047e-04
## Apobr 107.168314 -0.9809726 0.2410814 -4.069052 4.720479e-05
## Itgam 1551.719408 -0.8667869 0.2177043 -3.981486 6.848561e-05
## Itgax 76.749721 -1.7231101 0.3981678 -4.327597 1.507446e-05
## Dock1 67.448952 -1.0748446 0.3172980 -3.387492 7.053467e-04
## Ifitm6 663.805623 -1.2980036 0.2970194 -4.370096 1.241918e-05
## Jund 496.680233 0.9995620 0.1789544 5.585567 2.329384e-08
## Iqcn 70.123601 1.3592055 0.3291542 4.129389 3.637291e-05
## Ddx39a 238.392529 0.8363592 0.2040125 4.099549 4.139562e-05
## Ier2 339.428271 0.8108983 0.1760569 4.605888 4.107104e-06
## Ces1d 11.919910 -2.4537084 0.7239551 -3.389310 7.006868e-04
## St3gal2 88.117861 -0.9221118 0.2668888 -3.455041 5.502090e-04
## Crispld2 101.640268 -1.0926381 0.3063453 -3.566688 3.615218e-04
## Mmp12 37.434911 -5.3030321 1.2456682 -4.257179 2.070229e-05
## Mmp8 1160.436670 -0.8628719 0.1996189 -4.322595 1.542045e-05
## Sorl1 930.504301 -0.7350288 0.2169068 -3.388685 7.022858e-04
## Il10ra 116.034134 -0.9577016 0.2686355 -3.565060 3.637731e-04
## Pkm 970.124244 -0.7488659 0.1867988 -4.008943 6.099105e-05
## Pik3cb 135.444540 -1.0868643 0.2913901 -3.729929 1.915341e-04
## Ltf 7428.912233 -0.9584815 0.2789290 -3.436291 5.897363e-04
## Eomes 10.729335 -2.3358781 0.6775369 -3.447603 5.655851e-04
## Ccr2 333.437859 -1.5840455 0.2627479 -6.028766 1.652159e-09
## Tnfaip3 491.482503 1.6742184 0.2357063 7.102986 1.220900e-12
## Sgk1 501.865759 0.8730832 0.2353681 3.709437 2.077204e-04
## Adamts14 13.489327 -2.2325827 0.6130459 -3.641787 2.707516e-04
## Itgb2 1031.779954 -0.8807052 0.2351735 -3.744917 1.804536e-04
## Prss57 41.800254 -1.9456070 0.5716883 -3.403265 6.658564e-04
## Plxnc1 245.692246 -0.8926896 0.2098783 -4.253367 2.105798e-05
## Dusp6 357.384285 1.0230410 0.2016763 5.072687 3.922377e-07
## Lyz2 4421.351744 -1.0468202 0.2422598 -4.321064 1.552789e-05
## Ppm1h 212.789840 -0.8921082 0.2335941 -3.819052 1.339654e-04
## Ddit3 62.795830 1.1257726 0.2947776 3.819058 1.339621e-04
## Mgl2 29.903115 -1.6902948 0.4566262 -3.701703 2.141568e-04
## Bcl6b 36.944070 1.7286809 0.4507191 3.835384 1.253682e-04
## Slfn8 119.432808 -0.9136515 0.2533734 -3.605949 3.110148e-04
## Slfn4 530.056281 -0.8637450 0.2210971 -3.906633 9.359121e-05
## Ccl5 29.173569 -2.2995818 0.6063877 -3.792263 1.492807e-04
## Ccl4 25.625065 1.9472130 0.4656220 4.181961 2.890060e-05
## Wfdc17 49.245118 -1.8723226 0.4727177 -3.960762 7.471086e-05
## Arl5c 281.773850 0.8116742 0.2379413 3.411237 6.466894e-04
## Ace 57.141523 -1.3204677 0.3288271 -4.015690 5.927210e-05
## Klf11 59.277023 1.0326072 0.2976904 3.468729 5.229264e-04
## Gm36723 50.184361 -1.4814379 0.3852840 -3.845055 1.205255e-04
## Gm22513 821.522087 -1.1307338 0.2838874 -3.983036 6.804052e-05
## Zfp36l1 446.088731 0.9943620 0.2118156 4.694470 2.672987e-06
## Actn1 274.870254 -0.8699495 0.2125806 -4.092329 4.270627e-05
## Fos 978.275952 1.0956212 0.2647691 4.138026 3.503072e-05
## Galc 82.891061 -1.0295785 0.2952672 -3.486938 4.885848e-04
## Ighv8-7 10.643617 2.6476451 0.7793734 3.397146 6.809265e-04
## Gm35730 28.701648 1.9056577 0.5120409 3.721690 1.978935e-04
## Aoah 99.104914 -1.1526141 0.3254953 -3.541109 3.984497e-04
## Cmah 75.928543 -1.3456313 0.3759665 -3.579126 3.447456e-04
## Irf4 295.482623 0.9087470 0.2551446 3.561695 3.684690e-04
## Tgfbi 220.152432 -0.8824823 0.2518766 -3.503629 4.589636e-04
## Arrdc3 370.110199 0.8076989 0.1957495 4.126187 3.688279e-05
## Mef2c 1586.707689 1.1568150 0.3376243 3.426338 6.117787e-04
## Vcan 35.738525 -1.7084435 0.3705706 -4.610305 4.020786e-06
## Gm36161 62.754075 -1.4636888 0.3493906 -4.189262 2.798628e-05
## Vcl 482.922822 -0.7055499 0.2031617 -3.472849 5.149645e-04
## Bmpr1a 54.588182 -2.4221931 0.4980287 -4.863562 1.152921e-06
## Lgals3 361.613544 -1.4776667 0.3179289 -4.647790 3.355095e-06
## Mmp14 46.021690 -1.8139646 0.4534762 -4.000132 6.330714e-05
## Cebpe 193.839160 -1.2958437 0.3669742 -3.531158 4.137445e-04
## Ctsg 239.796481 -1.5894319 0.4117473 -3.860212 1.132888e-04
## Adra1a 9.178068 -5.8173314 1.2057644 -4.824600 1.402843e-06
## Dock5 285.151083 -0.8542187 0.2258004 -3.783070 1.549059e-04
## Fyb1 286.793529 -0.7696619 0.2225232 -3.458793 5.426010e-04
## Zfp622 167.776550 0.8253282 0.2426455 3.401374 6.704813e-04
## Ly6c2 436.171079 -1.2025387 0.3196160 -3.762449 1.682578e-04
## Il2rb 83.646898 -2.1444046 0.4521050 -4.743156 2.104138e-06
## Bin2 229.940714 -1.1062883 0.3063897 -3.610723 3.053443e-04
## Mefv 55.631663 -1.5746766 0.3831764 -4.109534 3.964585e-05
## Ciita 72.183446 -1.4018572 0.3918927 -3.577146 3.473665e-04
## Septin5 104.380152 -1.1784543 0.3349593 -3.518201 4.344829e-04
## Nfkbiz 427.098545 1.0680589 0.3037954 3.515717 4.385677e-04
## St3gal6 97.800821 0.9761665 0.2642685 3.693844 2.208895e-04
## Fpr2 170.810100 -1.0435345 0.2982143 -3.499278 4.665207e-04
## Rpl3l 32.655640 1.5290490 0.3851948 3.969548 7.200924e-05
## Prss34 165.843413 0.8643663 0.2192175 3.942962 8.048149e-05
## Dusp1 642.988058 0.8287506 0.2445521 3.388852 7.018600e-04
## H2-Aa 450.031532 -1.0939917 0.2797218 -3.910999 9.191526e-05
## H2-Eb1 466.341001 -1.1810818 0.2868497 -4.117424 3.831312e-05
## Pla2g7 116.415223 -1.2289662 0.3253515 -3.777349 1.585063e-04
## Treml2 236.649470 -0.6263603 0.1849694 -3.386291 7.084414e-04
## Tgif1 186.413028 0.8255509 0.2167275 3.809166 1.394363e-04
## Xdh 383.203840 -0.7620958 0.2039986 -3.735790 1.871268e-04
## Egr1 365.603106 1.7818740 0.3806285 4.681399 2.849231e-06
## Cd74 1305.077804 -1.1781185 0.2951437 -3.991677 6.560766e-05
## Ahnak 811.893401 -0.9999559 0.2767688 -3.612965 3.027153e-04
## Ms4a4c 42.239897 -1.7413902 0.3711365 -4.692048 2.704832e-06
## Ms4a4b 35.218239 -2.0047348 0.5721415 -3.503914 4.584726e-04
## Ms4a3 131.068001 -1.4609811 0.4209798 -3.470431 5.196242e-04
## Mpeg1 948.120915 -1.1322084 0.3059168 -3.701034 2.147224e-04
## Anxa1 1996.122060 -0.9565209 0.2339847 -4.087964 4.351761e-05
## Gda 653.009076 -0.9597831 0.2301804 -4.169700 3.050005e-05
## Myof 84.908920 -2.1194035 0.4046852 -5.237166 1.630613e-07
## Cybb 2218.810005 -1.0088288 0.2575414 -3.917152 8.960142e-05
## G6pdx 566.541587 -1.1233134 0.2599491 -4.321283 1.551247e-05
## Mir223hg 913.079303 -1.2106112 0.3000651 -4.034495 5.471981e-05
## Vsig4 12.265848 -2.7752944 0.7360879 -3.770330 1.630318e-04
## Rnf128 13.353695 -2.9442851 0.7892219 -3.730618 1.910108e-04
## Tmsb4x 2537.406467 -0.8949135 0.2543691 -3.518169 4.345350e-04
## Tlr8 87.760130 -1.0350131 0.2914918 -3.550745 3.841425e-04
## COX1 16567.633211 -0.8470856 0.2473184 -3.425081 6.146164e-04
## COX3 16.584291 -2.0395764 0.5112723 -3.989218 6.629152e-05
## padj
## Gm4956 1.897904e-02
## Gm41914 3.438463e-02
## Actr3 4.809468e-02
## Rab7b 3.889496e-02
## Fasl 9.834078e-03
## Xcl1 2.204708e-02
## Fcgr2b 3.889496e-02
## Fcgr4 4.992359e-02
## Pcp4l1 3.744312e-02
## Slamf7 4.905213e-02
## Atf3 1.346173e-02
## Olfm1 9.963647e-03
## Gm39861 2.381187e-03
## Rasgrp1 1.538223e-02
## Gatm 3.143157e-02
## A630026N12Rik 1.487362e-02
## Dusp2 1.342140e-02
## Tgif2 1.342140e-02
## Samhd1 9.834078e-03
## Pltp 3.821222e-02
## Zfp971 1.412426e-02
## Sirpb1a 3.769029e-02
## Zbtb7b 1.342140e-02
## S100a8 4.515612e-03
## Chil3 3.777903e-02
## S1pr1 3.881271e-02
## Cnn3 4.494075e-02
## Gem 2.627721e-02
## Atp6v0d2 3.021646e-02
## Klf4 1.538223e-02
## Slc31a2 3.881271e-02
## Megf9 4.687659e-02
## Jun 7.259311e-03
## Pde4b 4.160209e-03
## Zc3h12a 2.584813e-02
## Man1c1 4.498664e-02
## Arhgef10l 4.972193e-02
## Clnk 2.404181e-03
## Prom1 4.992359e-02
## Anxa3 2.949344e-02
## Hpse 4.596381e-03
## Btbd8 2.119962e-02
## Oasl2 4.905213e-02
## Ncf1 3.731204e-02
## Mgam 1.255888e-02
## Gimap3 2.075318e-02
## Gpnmb 6.990600e-04
## Skap2 4.992359e-02
## Gm44008 4.490987e-02
## Anxa4 4.972193e-02
## Rassf4 2.124832e-02
## C3ar1 4.494075e-02
## Clec4a3 2.446466e-02
## Klrb1a 4.239387e-02
## Cd69 1.336864e-03
## Clec7a 4.992359e-02
## Klre1 3.481885e-02
## Klrk1 1.336864e-03
## Klrc3 3.730017e-02
## Klri2 3.301522e-03
## Klra8 6.990600e-04
## Klra7 3.972616e-03
## Pirb 3.335989e-03
## Ncr1 2.137761e-03
## Pglyrp1 1.538223e-02
## Cd177 3.847578e-02
## Zfp36 1.336864e-03
## Rinl 2.443017e-02
## Ppp1r15a 2.252328e-03
## Svip 4.809468e-02
## Mfge8 1.136032e-02
## Anpep 2.820788e-02
## Gm15501 1.538223e-02
## Tenm4 2.527908e-02
## Bmal1 4.905213e-02
## Itpripl2 2.536392e-02
## Igsf6 1.930227e-02
## Apobr 1.136032e-02
## Itgam 1.343898e-02
## Itgax 5.825243e-03
## Dock1 4.992359e-02
## Ifitm6 5.280221e-03
## Jund 9.903762e-05
## Iqcn 1.022696e-02
## Ddx39a 1.077553e-02
## Ier2 2.381187e-03
## Ces1d 4.992359e-02
## St3gal2 4.498664e-02
## Crispld2 3.730017e-02
## Mmp12 7.259311e-03
## Mmp8 5.825243e-03
## Sorl1 4.992359e-02
## Il10ra 3.730017e-02
## Pkm 1.342140e-02
## Pik3cb 2.443017e-02
## Ltf 4.730872e-02
## Eomes 4.594929e-02
## Ccr2 1.053664e-05
## Tnfaip3 1.557257e-08
## Sgk1 2.536392e-02
## Adamts14 3.029330e-02
## Itgb2 2.397589e-02
## Prss57 4.992359e-02
## Plxnc1 7.259311e-03
## Dusp6 7.147130e-04
## Lyz2 5.825243e-03
## Ppm1h 2.010270e-02
## Ddit3 2.010270e-02
## Mgl2 2.559611e-02
## Bcl6b 1.930227e-02
## Slfn8 3.361859e-02
## Slfn4 1.538223e-02
## Ccl5 2.124832e-02
## Ccl4 9.451977e-03
## Wfdc17 1.401378e-02
## Arl5c 4.939235e-02
## Ace 1.342140e-02
## Klf11 4.446617e-02
## Gm36723 1.897904e-02
## Gm22513 1.343898e-02
## Zfp36l1 2.137761e-03
## Actn1 1.088367e-02
## Fos 1.015493e-02
## Galc 4.239387e-02
## Ighv8-7 4.992359e-02
## Gm35730 2.474639e-02
## Aoah 3.821222e-02
## Cmah 3.631688e-02
## Irf4 3.730017e-02
## Tgfbi 4.065334e-02
## Arrdc3 1.022696e-02
## Mef2c 4.809468e-02
## Vcan 2.381187e-03
## Gm36161 9.393817e-03
## Vcl 4.438090e-02
## Bmpr1a 1.336864e-03
## Lgals3 2.252328e-03
## Mmp14 1.342140e-02
## Cebpe 3.881271e-02
## Ctsg 1.829111e-02
## Adra1a 1.491106e-03
## Dock5 2.171235e-02
## Fyb1 4.494075e-02
## Zfp622 4.992359e-02
## Ly6c2 2.259083e-02
## Il2rb 2.064484e-03
## Bin2 3.328774e-02
## Mefv 1.053506e-02
## Ciita 3.631688e-02
## Septin5 3.930847e-02
## Nfkbiz 3.939388e-02
## St3gal6 2.584813e-02
## Fpr2 4.103774e-02
## Rpl3l 1.370863e-02
## Prss34 1.466488e-02
## Dusp1 4.992359e-02
## H2-Aa 1.538223e-02
## H2-Eb1 1.039753e-02
## Pla2g7 2.197552e-02
## Treml2 4.992359e-02
## Tgif1 2.068035e-02
## Xdh 2.443017e-02
## Egr1 2.137761e-03
## Cd74 1.342140e-02
## Ahnak 3.328563e-02
## Ms4a4c 2.137761e-03
## Ms4a4b 4.065334e-02
## Ms4a3 4.446617e-02
## Mpeg1 2.559611e-02
## Anxa1 1.088367e-02
## Gda 9.725703e-03
## Myof 5.199617e-04
## Cybb 1.538223e-02
## G6pdx 5.825243e-03
## Mir223hg 1.269002e-02
## Vsig4 2.212202e-02
## Rnf128 2.443017e-02
## Tmsb4x 3.930847e-02
## Tlr8 3.769029e-02
## COX1 4.809468e-02
## COX3 1.342140e-02Hier beginnen we met gebruiken van KEGG. KEGG ook wel Kyoto Encyclopedia of Genes and Genomes is een website die door middel van clusterprofiler gebruikt kan worden in Rstudio. clusterprofiler gebruikt de KEGG database voor info en “connect” met hun server om informatie te verkrijgen. Dit is te zien door een analyse te runnen en in de terminal is te zien dat er communicatie is tussen de server en jouw R.
Hier onder zoek ik de gene symbols uit. De gene symbols zijn nodig om entrez_id’s te vinden. De entrez id’s zijn nodig in latere stappen en deze entrez_id’s zijn te linken aan de gene symbolen.
#hier maak ik een data frame aan waar ik de gen namen in zet en de pvalues ook overneem
data_kegg <- data.frame(
#De gen naam is nodig voor het verkrijgen van de ENTREZID
Gene = row.names(deseq_shrink),
#Pvalue is nodig voor filtering in volgende stappen
pvalue = deseq_shrink$pvalue,
stringsAsFactors = FALSE
)
gene_symbols <- row.names(deseq_shrink)
head(gene_symbols)## [1] "Gm26206" "Xkr4" "Gm53491" "Rp1" "Sox17" "Gm22307"Met bitr kunnen we doormiddel van fromType en toType de gene symbolen omzetten naar de gewilde entrez_id’s verder is nog de Org.Mm.eg.db dit betekent Organism Mouse musculus eg database. We willen natuurlijk de entrez id’s van de muis terug krijgen en niet de mens daarom specificeren we de organism database. De entrez id’s zijn essentieel voor de KEGG analyse.
#hier gebruik ik bitr dit is een Biological Id TRanslator, hiermee vertaal ik mijn gen namen naar de ENTREID die de kegg tools verwachten met de data base van de mmu
entrez_ids <- bitr(rownames(resultaat_deseq_datafr), fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Mm.eg.db)## 'select()' returned 1:1 mapping between keys and columnshead(entrez_ids)## SYMBOL ENTREZID
## 1 Gm4956 241041
## 2 Gm41914 105246657
## 3 Actr3 74117
## 4 Rab7b 226421
## 5 Fasl 14103
## 6 Xcl1 16963#entrez_ids$ENTREID geeft aan de gene de juiste id's mee voor elk gen, orgbd vraagt welke database org.Mm.eg.db is de database van de mus musculus, keytipe geeft aan dat de genen die worden gebruikt in de analyse geïdentificeerd met entrez ID's. ont = BP geeft aan dat wel Biologische processen willen. Door dit te veranderen kan je ook MF - moleculaire functie krijgen en CC - cellular component. pAdjustMethod geeft aan dat de Benjamini-Hochberg (BH) methode wordt gebruikt om de p-waarden voor meerdere tests te corrigeren. Deze methode helpt om het risico op vals-positieve resultaten te verminderen. qvalueCutoff = 0.05 stelt dat alles onder 0.05 statistisch significant is en dus belangrijk is. readable = true zet de gen symbolen inplaats van de entrez IDs zodat het te bekijken is.
go_results <- enrichGO(gene = entrez_ids$ENTREZID,
OrgDb = org.Mm.eg.db,
keyType = "ENTREZID",
ont = "BP",
pAdjustMethod = "BH",
qvalueCutoff = 0.05,
readable = TRUE)
head(go_results)## ID Description GeneRatio BgRatio
## GO:0002697 GO:0002697 regulation of immune effector process 25/168 473/28905
## GO:0045088 GO:0045088 regulation of innate immune response 24/168 486/28905
## GO:0097529 GO:0097529 myeloid leukocyte migration 18/168 255/28905
## GO:0007159 GO:0007159 leukocyte cell-cell adhesion 22/168 445/28905
## GO:0002274 GO:0002274 myeloid leukocyte activation 18/168 281/28905
## GO:0030595 GO:0030595 leukocyte chemotaxis 17/168 242/28905
## RichFactor FoldEnrichment zScore pvalue p.adjust
## GO:0002697 0.05285412 9.093741 13.57029 5.797906e-17 1.784595e-13
## GO:0045088 0.04938272 8.496473 12.74335 1.165385e-15 1.793528e-12
## GO:0097529 0.07058824 12.144958 13.66782 1.219735e-14 1.251448e-11
## GO:0007159 0.04943820 8.506019 12.20071 1.921520e-14 1.478610e-11
## GO:0002274 0.06405694 11.021225 12.90689 6.524299e-14 3.816602e-11
## GO:0030595 0.07024793 12.086408 13.24191 7.439769e-14 3.816602e-11
## qvalue
## GO:0002697 1.171177e-13
## GO:0045088 1.177039e-12
## GO:0097529 8.212885e-12
## GO:0007159 9.703677e-12
## GO:0002274 2.504722e-11
## GO:0030595 2.504722e-11
## geneID
## GO:0002697 Xcl1/Fcgr2b/Rasgrp1/Zbtb7b/Zc3h12a/Clnk/Ncf1/Gimap3/Klrb1a/Clec7a/Klre1/Klrk1/Klrc3/Klri2/Pglyrp1/Cd177/Itgam/Ccr2/Tnfaip3/Itgb2/Lgals3/Nfkbiz/Cd74/Anxa1/Vsig4
## GO:0045088 Rab7b/Xcl1/Rasgrp1/Samhd1/S100a8/Clnk/Ncf1/Gimap3/Klrb1a/Clec7a/Klre1/Klrk1/Klrc3/Klri2/Mmp12/Ltf/Tnfaip3/Irf4/Mefv/Nfkbiz/Fpr2/Cd74/Vsig4/Tlr8
## GO:0097529 Xcl1/S100a8/Pde4b/C3ar1/Cd177/Itgam/Ccr2/Itgb2/Ccl5/Ccl4/Lgals3/Mmp14/Ctsg/Fpr2/Dusp1/Pla2g7/Cd74/Anxa1
## GO:0007159 Xcl1/Rasgrp1/Sirpb1a/Zbtb7b/S100a8/Klf4/Zc3h12a/Gimap3/Gpnmb/Cd177/Itgam/Ccr2/Itgb2/Ccl5/Lgals3/Ctsg/Nfkbiz/H2-Aa/H2-Eb1/Cd74/Anxa1/Vsig4
## GO:0002274 Fcgr4/Rasgrp1/Jun/Clnk/Anxa3/Gimap3/Klrk1/Pirb/Cd177/Itgam/Jund/Mmp8/Ccr2/Itgb2/Ccl5/Irf4/Ctsg/Anxa1
## GO:0030595 Xcl1/S100a8/S1pr1/Pde4b/C3ar1/Itgam/Ccr2/Itgb2/Ccl5/Ccl4/Lgals3/Ctsg/Fpr2/Dusp1/Pla2g7/Cd74/Anxa1
## Count
## GO:0002697 25
## GO:0045088 24
## GO:0097529 18
## GO:0007159 22
## GO:0002274 18
## GO:0030595 17De GO-verrijkingsanalyse helpt te begrijpen welke biologische processen, functies of cellulaire componenten oververtegenwoordigd zijn in een bepaalde genenlijst. Dit geeft inzicht in biologische context van de genen die bestuurt zijn
Het uitvoeren van een GO-verrijkingsanalyse helpt onderzoekers te begrijpen welke biologische processen, functies of cellulaire componenten oververtegenwoordigd zijn in een bepaalde genenlijst. Dit kan waardevolle inzichten opleveren over de biologische context van de genen die je hebt bestudeerd.
De volgende data wordt terug gegeven GeneRatio - het percentage van het aantal genen dat aanwezig is in deze GO-term ten opzichte van het totale aantal genen in deze categorie.
BgRatio - verhouding van alle genen die in deze term zijn geannoteerd
RichFactor - De verhouding vanhet differentieel tot expressie gebrachte eiwitnummer dat in deze padterm is geannoteerd tot alle eiwitnummers die zijn geannoteerd
FoldEnrichment - het percentage genen in uw lijst dat tot een pad behoort, gedeeld door het overeenkomstige percentage in de achtergrond
zScore - geeft het aantal standaarddeviaties weer dat een gegevenspunt van het gemiddelde verwijderd is als een standaardscore ID - geeft het Gene Ontology ID weer wat kan worden opgezocht om te zien wat dit gen doet - bijvoorbeeld GO:0051051 in google geeft een website met uitleg wat het doet voor functie : negative regulation of transport https://www.informatics.jax.org/vocab/gene_ontology/GO:0051051
Description - geeft ook weer wat er te vinden is over de gen functie
geneID - geeft welke genen bij deze functie werken.
#entrez_ids$ENTREID geeft aan de gene de juiste id's mee voor elk gen, orgbd vraagt welke database org.Mm.eg.db is de database van de mus musculus,
ego <- enrichGO(gene = entrez_ids$ENTREZID,
universe = names(entrez_ids$SYMBOL),
OrgDb = org.Mm.eg.db,
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
head(ego)## ID Description GeneRatio BgRatio
## GO:0002102 GO:0002102 podosome 4/170 30/28618
## GO:0070062 GO:0070062 extracellular exosome 6/170 106/28618
## GO:0042613 GO:0042613 MHC class II protein complex 3/170 13/28618
## GO:0045121 GO:0045121 membrane raft 9/170 305/28618
## GO:0098857 GO:0098857 membrane microdomain 9/170 306/28618
## GO:1903561 GO:1903561 extracellular vesicle 6/170 128/28618
## RichFactor FoldEnrichment zScore pvalue p.adjust
## GO:0002102 0.13333333 22.445490 9.084800 2.919493e-05 0.004417486
## GO:0070062 0.05660377 9.528746 6.800407 4.199625e-05 0.004417486
## GO:0042613 0.23076923 38.847964 10.551249 5.637809e-05 0.004417486
## GO:0045121 0.02950820 4.967445 5.384910 9.291897e-05 0.004417486
## GO:0098857 0.02941176 4.951211 5.371756 9.525372e-05 0.004417486
## GO:1903561 0.04687500 7.890993 6.040198 1.199213e-04 0.004417486
## qvalue geneID Count
## GO:0002102 0.003368494 Actr3/Vcl/Dock5/Bin2 4
## GO:0070062 0.003368494 Fasl/Prom1/Anpep/Ace/Ahnak/Anxa1 6
## GO:0042613 0.003368494 H2-Aa/H2-Eb1/Cd74 3
## GO:0045121 0.003368494 Fasl/Hpse/Cd177/Itgam/Itgb2/Vcl/Bmpr1a/Ahnak/Myof 9
## GO:0098857 0.003368494 Fasl/Hpse/Cd177/Itgam/Itgb2/Vcl/Bmpr1a/Ahnak/Myof 9
## GO:1903561 0.003368494 Fasl/Prom1/Anpep/Ace/Ahnak/Anxa1 6Hier word de enrichKEGG gemaakt. De enrich kegg zoekt uit welke genen statistich significant zijn in hun expressie, er zijn 3 dingen nodig om de analyse te doen. Er moeten entrez_ids zijn deze zijn gemaakt in een eerdere functie, organisme moet worden gespecificeerd en de pvalueCutoff moet worden aangegeven die weer 0.05 is.
De resultaten zet ik in de variabele kk en deze zet ik om in een dataframe zodat deze makkelijker te lezen is.
kk <- enrichKEGG(gene = entrez_ids$ENTREZID,
organism = 'mmu',
pvalueCutoff = 0.05)## Reading KEGG annotation online: "https://rest.kegg.jp/link/mmu/pathway"...## Reading KEGG annotation online: "https://rest.kegg.jp/list/pathway/mmu"...kk## #
## # over-representation test
## #
## #...@organism mmu
## #...@ontology KEGG
## #...@keytype kegg
## #...@gene chr [1:181] "241041" "105246657" "74117" "226421" "14103" "16963" "14130" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...38 enriched terms found
## 'data.frame': 38 obs. of 14 variables:
## $ category : chr "Human Diseases" "Human Diseases" "Cellular Processes" "Organismal Systems" ...
## $ subcategory : chr "Infectious disease: parasitic" "Infectious disease: bacterial" "Transport and catabolism" "Development and regeneration" ...
## $ ID : chr "mmu05140" "mmu05152" "mmu04145" "mmu04380" ...
## $ Description : chr "Leishmaniasis - Mus musculus (house mouse)" "Tuberculosis - Mus musculus (house mouse)" "Phagosome - Mus musculus (house mouse)" "Osteoclast differentiation - Mus musculus (house mouse)" ...
## $ GeneRatio : chr "9/103" "12/103" "11/103" "9/103" ...
## $ BgRatio : chr "70/9793" "180/9793" "183/9793" "136/9793" ...
## $ RichFactor : num 0.1286 0.0667 0.0601 0.0662 0.0805 ...
## $ FoldEnrichment: num 12.22 6.34 5.72 6.29 7.65 ...
## $ zScore : num 9.72 7.45 6.64 6.41 6.42 ...
## $ pvalue : num 4.22e-08 3.68e-07 3.22e-06 1.23e-05 3.40e-05 ...
## $ p.adjust : num 9.04e-06 3.93e-05 2.30e-04 6.56e-04 1.46e-03 ...
## $ qvalue : num 6.80e-06 2.96e-05 1.73e-04 4.94e-04 1.10e-03 ...
## $ geneID : chr "246256/16476/17969/16409/16414/14281/14960/14969/13058" "14130/246256/242341/56644/16409/16411/16154/16414/12265/14960/14969/16149" "226421/14130/246256/242341/17969/56644/16409/16414/14960/14969/13058" "14130/246256/320832/16476/17969/18733/16478/74769/14281" ...
## $ Count : int 9 12 11 9 7 8 8 10 8 7 ...
## #...Citation
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
## clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
## The Innovation. 2021, 2(3):100141kk_data_frame <- as.data.frame(kk)hier is de module over-representation analysis kegg analyse te vinden. Met de pvalueCutoff en de qvalueCutoff
Er zijn verschillende dingen te zien in deze dataframe
GeneRatio - het percentage van het aantal genen dat aanwezig is in deze GO-term ten opzichte van het totale aantal genen in deze categorie.
Bgratio - verhouding van alle genen die in deze term zijn geannoteerd
RichFactor - De verhouding vanhet differentieel tot expressie gebrachte eiwitnummer dat in deze padterm is geannoteerd tot alle eiwitnummers die zijn geannoteerd
FoldEnrichment - het percentage genen in uw lijst dat tot een pad behoort, gedeeld door het overeenkomstige percentage in de achtergrond
zScore - geeft aan hoeveel standaarddeviaties een score van het gemiddelde af zit
pvalue - is een getal tussen 0 en 1, waarmee je bepaalt of een steekproefuitkomst statistisch significant is.
geneID - geeft welke genen bij deze functie werken.
count - hoeveel er aanwezig zijn die voldoen aan de paramaters
mkk <- enrichMKEGG(gene = entrez_ids$ENTREZID,
organism = 'mmu',
pvalueCutoff = 0.05,
qvalueCutoff = 0.05)## Reading KEGG annotation online: "https://rest.kegg.jp/link/mmu/module"...## Reading KEGG annotation online: "https://rest.kegg.jp/list/module"...mkk_data_frame <- as.data.frame(mkk)
mkk_data_frame## ID Description GeneRatio
## M00959 M00959 Guanine ribonucleotide degradation, GMP => Urate 2/10
## BgRatio RichFactor FoldEnrichment zScore pvalue p.adjust
## M00959 12/886 0.1666667 14.76667 5.127236 0.007130478 0.04991335
## qvalue geneID Count
## M00959 0.03752883 22436/14544 2Er is te zien dat er maar 1 gen uit de enrichMKEGG komt en dit is de M00959, deze gen heeft te maken met de Huanine ribonucleotide degredation. Deze gen is niet vernoemd in ons onderzoek waar wij ons onderzoek op baseren en om deze reden zullen wij niks verder doen met de data die hieruit is gekomen.
23-10-2024
#hier zet ik een lege lijst in de kk_list_id
kk_list_id <- list()
#hier neem ik alle ID's die uit de kegg analyse komen die belangrijk zijn en zet deze in de nieuwe lijst
kk_list_id <- kk_data_frame$ID
kk_list_id## [1] "mmu05140" "mmu05152" "mmu04145" "mmu04380" "mmu05323" "mmu04650"
## [7] "mmu05150" "mmu04613" "mmu05135" "mmu05146" "mmu04928" "mmu04670"
## [13] "mmu04666" "mmu04061" "mmu04810" "mmu05142" "mmu04620" "mmu04659"
## [19] "mmu05162" "mmu05418" "mmu05166" "mmu04010" "mmu04668" "mmu04932"
## [25] "mmu04148" "mmu04662" "mmu05417" "mmu04658" "mmu05202" "mmu04657"
## [31] "mmu04610" "mmu05221" "mmu04614" "mmu05415" "mmu05020" "mmu05171"
## [37] "mmu05100" "mmu05133"De reden om deze lijst te maken is het makkelijk bekijken van welke gen pathways interresant zijn.
De uiteindlijk interresante in mijn analyse zijn de volgende pathways.
#mmu04666: B-celreceptoren
#mmu04612: B-celactivatie
#mmu04611: B-celdifferentiatie
#mmu04670: B-celontwikkeling
#hier word de bijbehorende pathway geopend in de KEGG website.
browseKEGG(kk, "mmu04611")Hier worden de pathway’s met de expressie data aangeroepen wat hier uitkomt zijn foto’s. Deze foto zoals genoemd in de inleiding hebben een bepaalde structuur.
het belangerijkste om te noemen is de expressie schaal wat dus de kleuren schaal is. De kleuren schaal en cijfers betekenen ook wat, wanneer de gen een witte/grijze kleur heeft is de verandering in de expressie 0 wat dus betekent dat de expressie van dit gen normaal is. Wanneer een gen block een groene kleur heeft en dus richting de -1 gaat zegt dit dat de expressie is minder is dan de norm. Wanneer een gen een rode kleur heeft en dus richting de 1 gaat zegt dit dat het een overexpressie heeft meer dan het normaal.
#hierw worden alle pathviews
library(pathview)## ## ##############################################################################
## Pathview is an open source software package distributed under GNU General
## Public License version 3 (GPLv3). Details of GPLv3 is available at
## http://www.gnu.org/licenses/gpl-3.0.html. Particullary, users are required to
## formally cite the original Pathview paper (not just mention it) in publications
## or products. For details, do citation("pathview") within R.
##
## The pathview downloads and uses KEGG data. Non-academic uses may require a KEGG
## license agreement (details at http://www.kegg.jp/kegg/legal.html).
## ##############################################################################mmu04612_pathview <- pathview(gene.data = resultaat_deseq_datafr,
pathway.id = "mmu04612",
species = "mmu")## Warning: None of the genes or compounds mapped to the pathway!
## Argument gene.idtype or cpd.idtype may be wrong.## Warning in colnames(plot.data)[c(1, 3, 9:ncs)] <- c("kegg.names", "all.mapped",
## : number of items to replace is not a multiple of replacement length## 'select()' returned 1:1 mapping between keys and columns## Info: Working in directory /Users/jarnoduiker/github_bioinf/Genomics_transcriptomics_analysis/analysis/Jarno## Info: Writing image file mmu04612.pathview.pngmmu04666_pathview <- pathview(gene.data = resultaat_deseq_datafr,
pathway.id = "mmu04666",
species = "mmu")## Warning: None of the genes or compounds mapped to the pathway!
## Argument gene.idtype or cpd.idtype may be wrong.## Warning in colnames(plot.data)[c(1, 3, 9:ncs)] <- c("kegg.names", "all.mapped",
## : number of items to replace is not a multiple of replacement length## 'select()' returned 1:1 mapping between keys and columns## Info: Working in directory /Users/jarnoduiker/github_bioinf/Genomics_transcriptomics_analysis/analysis/Jarno## Info: Writing image file mmu04666.pathview.pngmmu04611_pathview <- pathview(gene.data = resultaat_deseq_datafr,
pathway.id = "mmu04611",
species = "mmu")## Warning: None of the genes or compounds mapped to the pathway!
## Argument gene.idtype or cpd.idtype may be wrong.## Warning in colnames(plot.data)[c(1, 3, 9:ncs)] <- c("kegg.names", "all.mapped",
## : number of items to replace is not a multiple of replacement length## 'select()' returned 1:1 mapping between keys and columns## Info: Working in directory /Users/jarnoduiker/github_bioinf/Genomics_transcriptomics_analysis/analysis/Jarno## Info: Writing image file mmu04611.pathview.pngmmu04670_pathview <- pathview(gene.data = resultaat_deseq_datafr,
pathway.id = "mmu04670",
species = "mmu")## Warning: None of the genes or compounds mapped to the pathway!
## Argument gene.idtype or cpd.idtype may be wrong.## Warning in colnames(plot.data)[c(1, 3, 9:ncs)] <- c("kegg.names", "all.mapped",
## : number of items to replace is not a multiple of replacement length## 'select()' returned 1:1 mapping between keys and columns## Info: Working in directory /Users/jarnoduiker/github_bioinf/Genomics_transcriptomics_analysis/analysis/Jarno## Info: Writing image file mmu04670.pathview.pngDe visualisaties zijn te vinden in de github map, deze hebben dus incorrecte schaal en de expressie hierin is dus niet goed.
Floris zijn versie van de KEGG pathway - problemen met visualisatie opgelost.
Deseq - floris
Statistische bewerkingen die door DESeq2 uitgevoerd
worden:
Een van de tests zorgt er voor dat er rekening gehouden wordt met de false discovery rate (FDR). Hiervoor is de Benjamini Hockberg compensatie uitgevoerd om hier rekening mee te houden. De false discovery rate is het aantal genen in dit geval die ondanks dat ze niet echt differentialy expressed zijn toch zo lijken. Dit zijn dus false positieven. Wanneer de p-waarde bijvoorbeeld 0.01 is en er 4000 genen zijn, zijn er 40 false positieven bij. De kans is namelijk 1% met een p-waarde van 0.01 of lager.
levels(coldata$Groep) geeft “Jong” als eerste factor, dat betekend dat up- en downregulatie tenopzichte van de jonge groep bepaald zal worden.
dds <- DESeqDataSetFromMatrix(countData = counts_df,
colData = coldata,
design= ~ Groep)
dds <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingresultsNames(dds) # De naam is: Groep_Oud_vs_Jong## [1] "Intercept" "Groep_Oud_vs_Jong"res <- results(dds, name="Groep_Oud_vs_Jong")
# Gebruikt lfcschrink om de data met Log2 te transformeren:
res <- lfcShrink(dds, coef="Groep_Oud_vs_Jong", type="apeglm")## using 'apeglm' for LFC shrinkage. If used in published research, please cite:
## Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
## sequence count data: removing the noise and preserving large differences.
## Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895summary(res)##
## out of 28500 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up) : 1587, 5.6%
## LFC < 0 (down) : 1142, 4%
## outliers [1] : 3, 0.011%
## low counts [2] : 12953, 45%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?resultsDe betekenis van de kolommen van het resultaat van DESeq2:
mcols(res)$description## [1] "mean of normalized counts for all samples"
## [2] "log2 fold change (MAP): Groep Oud vs Jong"
## [3] "posterior SD: Groep Oud vs Jong"
## [4] "Wald test p-value: Groep Oud vs Jong"
## [5] "BH adjusted p-values"Visualisatie van Differentially expressed genes:
MA-Plot
plotMA(res, ylim=c(-2,2))
Figuur 18: Ma plot 2de deseq
Hier is een MA-plot te zien waar dus nog heel veel ruis in zit waardoor er geen goede conclusies uit zijn te halen.
plotCounts(dds, gene=which.min(res$padj), intgroup="Groep")
Figuur 19: gen met kleinste p waarde
Visualisatie van het gen met de kleinste p-waarde: dat is: Gpnmb. - er is hier te zien dat oude muizen dit gen dus hoog expressen en de jonge muizen amper.
hieronder is het gen te zien waarop deze grafiek is gebaseerd.
counts_df[which.min(res$padj),]## SRR21754408 SRR21754417 SRR21754418 SRR21754419 SRR21754420 SRR21754421
## Gpnmb 25 50 1 6 277 816
## SRR21754422 SRR21754423
## Gpnmb 1008 697Plot van relevante genen:
Efb1, pax5 en foxo5
d <- plotCounts(dds, gene="Ebf1", intgroup="Groep",
returnData=TRUE)
ggplot(d, aes(x=Groep, y=count)) +
geom_point(position=position_jitter(w=0.1,h=0)) +
scale_y_log10(breaks=c(25,100,400))
Figuur 20: Ebf1 counts nieuwe deseq
Hier is weer een genexpressie plot te zien van gen Ebf1. er is weer te zien dat de jonge muizen hoog in de expressie en de oude laag maar de uitschieter is weg bij de oude muizen.
d <- plotCounts(dds, gene="Pax5", intgroup="Groep",
returnData=TRUE)
ggplot(d, aes(x=Groep, y=count)) +
geom_point(position=position_jitter(w=0.1,h=0)) +
scale_y_log10(breaks=c(25,100,400))
Figuur 21: Pax5 count nieuwe deseq
Hier is weer een genexpressie plot te zien van gen Pax5. er is weer te zien dat de jonge muizen hoog in de expressie en de oude laag maar de uitschieter is weg bij de oude muizen. verder is er een oude muis die heel lage expressie in het algemeen heeft dus deze vergelijken met de andere samples zou ook nog interresant kunnen zijn.
# Plot gevonden op de website: https://introtogenomics.readthedocs.io/en/latest/2021.11.11.DeseqTutorial.html
vsd <- vst(dds)
#head(assay(vsd), 30)
select <- order(rowMeans(counts(dds,normalized=TRUE)),
decreasing=TRUE)[1:10]
df <- as.data.frame(colData(dds))
pheatmap(assay(vsd)[select,], annotation_col=df)
Figuur 22: heatmap genen
Dit is een pheatmap plot. Het is een visualisatie van gegevens in de vorm van een hittekaart, waarbij de waarden van een matrix visueel worden weergegeven met kleuren. Dit is handig om patronen, correlaties en trends in data te identificeren, vooral wanneer je te maken hebt met veel variabelen of een grote dataset.
Uit deze heatmap is te halen dat jong gesplit is bij de gm26917 gen dit is te zien door de donkere kleur bij de ene helft jonge muozen en de andere jonge muizen zijn geel. verder zou hier nog uit te halen zijn dat er 2 jonge muizen zijn met hoge expressie voor deze genen en de andere 2 jonge muizen lage gen epxressie hebben voor deze genen.
sampledists <- dist(t(assay(vsd)))
# Convert the 'dist' object into a matrix for creating a heatmap
sampleDistMatrix <- as.matrix(sampledists)
pheatmap(sampleDistMatrix, show_colnames = FALSE,
annotation_col = coldata, # Gebruikt dezelfde coldata die ook voor DESeq2 gebruikt is.
clustering_distance_rows = sampledists,
clustering_distance_cols = sampledists,
main = "Euclidean Sample Distances")
Figuur 23: Euclidean sample distance
Euclidean Sample Distances is een manier om de afstand te berekenen tussen twee punten in een ruimte. In de context van sample distances (afstanden tussen monsters of datapunten) verwijst het naar de afstand tussen twee gegevenspunten.
er is te zien dat de sample 4418 en 4420 een donker rode vlak hebben wat kan beteken dat ze ver uit elkaar liggen.
Principal component analysis:
Hier is een pcaplot te zien, een PCA plot is een grafische weergave van de resultaten van Principal Component Analysis (PCA), een techniek die wordt gebruikt om de meest significante patronen in een dataset te identificeren. De plot toont meestal de eerste twee of drie principale componenten als assen, waardoor we de data kunnen visualiseren in lagere dimensies. Dit helpt om structuren, clusters of trends in de data te ontdekken, zelfs wanneer de oorspronkelijke dataset veel variabelen bevat.
plotPCA(vsd, intgroup=c("Groep"))## using ntop=500 top features by variance
Figuur 24: PCA
In onze pcaplot is te zien dat er enorm grote variatie zit tussen de jonge samples, wel 15-20 punten.
degPlot(dds = dds, res = res, n = 6, xs = "Groep")## No genes were mapped to rowData. check ann parameter values.## Using gene as id variables## `geom_smooth()` using formula = 'y ~ x'
Figuur 25: top expressed genen
De degPlot() function weergeeft volgens de documentatie enkele “top genen” en de expressie verschillen tussen jong en oud. opmerkelijk is dat alle genen in de oude muizen hogere expressie hebben.
deseq.volcano <- function(res, datasetName) {
return(EnhancedVolcano(res, x = 'log2FoldChange', y = 'padj',
lab=rownames(res),
title = paste(datasetName, "Jonge vs Oude muizen"),
subtitle = bquote(italic('FDR <= 0.05 and absolute FC >= 2')),
# Change text and icon sizes
labSize = 3, pointSize = 1.5, axisLabSize=10, titleLabSize=12,
subtitleLabSize=8, captionLabSize=10,
# Disable legend
legendPosition = "none",
# Set cutoffs
pCutoff = 0.05, FCcutoff = 2))
}
deseq.volcano(res = res, datasetName = "Volcano plot:")
Figuur 26: Volcano plot: Jonge vs Oude muizen
De analyse van floris
Een volcano plot weergeeft alle genen maar enkel de genen die statistisch significant zijn en een hoge fold change hebben worden rood gekleurd. De genen die wel statistisch significant zijn maar geen significante fold change hebben zijn blauw gekleurd. Daarnaast zijn alle genen die geen statistisch significante p-waarde hebben en ook geen significante fold change hebben grijs gekleurd. De -log10(p-waarde) is op de y-as weergeven en op de x-as staat de log2(fold change).
Eigen analyse
Er zijn dus een paar genen te zien die wel interessant zouden zijn om te onderzoeken wat voor impact ze hebben op deze muizen. dit zijn bijvoorbeeld Mmp12, Hbb-bs of Fosb. Verder zijn de pax5 en ebf1 genen niet te zien in deze volcano plot.
Pathview:
Het probleem bij de eerste pathview is dat de expressies niveaus op een juiste manier hier in krijgen niet zo goed ging. Hierom heeft floris het gehele process nogmaals doorgelopen om te zien of hij de bugs kon fixen en het gehele process efficienter maken. Dit is geluk en er zijn een paar veranderingen. Floris heeft een pval_threshold van 0.05 hierdoor komen alleen statistisch significante resultaten in de deseq.degs.logfc
Van een aantal relevante pathways zijn afbeeldingen gemaakt, de Kegg pathway codes zijn afkomstig van de volgende website: https://www.genome.jp/kegg-bin/show_organism?menu_type=pathway_maps&org=mmu
De data die gebruikt wordt hier is de data van de deseq resultaten.
pval_threshold <- 0.05
data(gene.idtype.list);
head(gene.idtype.list)## [1] "SYMBOL" "GENENAME" "ENSEMBL" "ENSEMBLPROT" "UNIGENE"
## [6] "UNIPROT"#hier word een subset gemaakt met de logfoldchange met 1 voorwaarde dat de padj waarde onder de pval_threshold ligt.
deseq.degs.logfc <- subset(res, padj < pval_threshold, select = log2FoldChange)Het verschil met de pathway die ik gemaakt heb en die van floris is dus dat floris de logfc subset heeft en deze gebruikt als gene.data, in mijn pathview is de subset gebruikt
hier is nog een pathview voorbeeld van mij.
mmu04670_pathview <- pathview(gene.data = resultaat_deseq_datafr, pathway.id = “mmu04670”, species = “mmu”)
Als we deze vergelijken met die van floris is te zien dat ik de gehele subset mee geef en dat floris voor gene.data alleen een subset meegeeft met de logfoldchange waardoor de juiste expressie levels dus te zien zijn (-1 tot 1) ook is te zien dat er cpd.data= res toegevoegd is. De cpd staat voor compound en deze word terug gelinked aan de resultaten van de deseq.
Deze twee verandering hebben gezorgt dat de pathview er uit ziet zoals het moet en daarom zijn deze pathview resultaten leidend en gebruiken we deze voor de resultaten.
T cell receptor signaling pathway:
# 04660 T cell receptor signaling pathway
pathview(gene.data=deseq.degs.logfc,
cpd.data=res,
,gene.idtype = "SYMBOL",
pathway.id="04660",
species="mmu")## 'select()' returned 1:1 mapping between keys and columns## [1] "Note: 1 of 2119 unique input IDs unmapped."## 'select()' returned 1:1 mapping between keys and columns## Warning: None of the genes or compounds mapped to the pathway!
## Argument gene.idtype or cpd.idtype may be wrong.## Warning in colnames(plot.data)[c(1, 3, 9:ncs)] <- c("kegg.names", "all.mapped",
## : number of items to replace is not a multiple of replacement length## Info: Working in directory /Users/jarnoduiker/github_bioinf/Genomics_transcriptomics_analysis/analysis/Jarno## Info: Writing image file mmu04660.pathview.png# Om de afbeelding te weergeven in R:
knitr::include_graphics("mmu04660.pathview.png")
Figuur 27: T cell receptor signaling pathway
T cell receptor signaling pathway: Er is overlap tussen de surface eiwitten van B- en T-lymfocyten. Dus ondanks dat de gesequencede cellen B-cellen zijn en geen T-cellen is het volgende pathway ook interesant. hier kunnen met over en onderexpressie van essentiele genen dus ook problemen ontstaan bij dingen zoals imuumsysteem reactie
# 04662 B cell receptor signaling pathway
pathview(gene.data=deseq.degs.logfc,
cpd.data=res,
,gene.idtype = "SYMBOL",
pathway.id="04662",
species="mmu")## 'select()' returned 1:1 mapping between keys and columns## [1] "Note: 1 of 2119 unique input IDs unmapped."## 'select()' returned 1:1 mapping between keys and columns## Warning: None of the genes or compounds mapped to the pathway!
## Argument gene.idtype or cpd.idtype may be wrong.## Warning in colnames(plot.data)[c(1, 3, 9:ncs)] <- c("kegg.names", "all.mapped",
## : number of items to replace is not a multiple of replacement length## Info: Working in directory /Users/jarnoduiker/github_bioinf/Genomics_transcriptomics_analysis/analysis/Jarno## Info: Writing image file mmu04662.pathview.png# Om de afbeelding te weergeven in R:
knitr::include_graphics("mmu04662.pathview.png")
Figuur 28: B cell receptor signaling pathway
B cell receptor signaling pathway: De Genen die belangrijk zijn in ons originele onderzoek hebben afkomst van de B-cellen waardoor de B-cell receptoren pathway extreem relevant is. Wanneer er grote expressie verschillen zijn tussen de jonge en oude groep en er meer voorloper-b-cellen aanwezig zijn bij 1 van de groepen kan dit verklaard worden door expressie verschillen van de genen die verantwoordelijk zijn voor het maken van receptoreiwitten die als markers worden gebruikt. Hierdoor kunnen we ook zien dat wanneer er onder of overexpressie.
in onze pathway is te zien dat er een paar genen over expresses zijn en andere underexpressed. dit zou dus impact kunnen hebben op hoe sommige signalen worden verstuurd. Dit zou voor verschillende problemen kunnen zorgen. zoals problemen in imuum reactie.
# 04211
pathview(gene.data=deseq.degs.logfc,
cpd.data=res,
,gene.idtype = "SYMBOL",
pathway.id="04211",
species="mmu")## 'select()' returned 1:1 mapping between keys and columns## [1] "Note: 1 of 2119 unique input IDs unmapped."## 'select()' returned 1:1 mapping between keys and columns## Warning: None of the genes or compounds mapped to the pathway!
## Argument gene.idtype or cpd.idtype may be wrong.## Warning in colnames(plot.data)[c(1, 3, 9:ncs)] <- c("kegg.names", "all.mapped",
## : number of items to replace is not a multiple of replacement length## Info: Working in directory /Users/jarnoduiker/github_bioinf/Genomics_transcriptomics_analysis/analysis/Jarno## Info: Writing image file mmu04211.pathview.png# Om de afbeelding te weergeven in R:
knitr::include_graphics("mmu04211.pathview.png")
Figuur 29: Longevity regulating pathway
Analyse van floris
Longevity regulating pathway: De twee groepen die vergeleken worden zijn jonge en oude muizen, daarom is het longevity pathway dat enkele genen weergeeft die betrokken zijn bij het verouderingsproces ook belangrijk. Niet alle genen die betrokken zijn bij dit proces zijn hierop afgebeeld omdat het niet bekend is welke genen verantwoordelijk zijn. Foxo is downregulated in de oude muizengroep tegenover de jonge groep. Van dit gen is bekend dat het levensduur gerelateerd is. Bron: https://cellandbioscience.biomedcentral.com/articles/10.1186/s13578-021-00700-7
Analyse van mij
Problemen in deze pathway zullen dus impact hebben op de levensduur van de in de geval muis. Al is dit niet met zekerheid te zeggen omdat het nog niet bekend is welke genen specifiek belangrijk zijn voor levensduur.
# 04640
pathview(gene.data=deseq.degs.logfc,
cpd.data=res,
,gene.idtype = "SYMBOL",
pathway.id="04640",
species="mmu")## 'select()' returned 1:1 mapping between keys and columns## [1] "Note: 1 of 2119 unique input IDs unmapped."## 'select()' returned 1:1 mapping between keys and columns## Info: Working in directory /Users/jarnoduiker/github_bioinf/Genomics_transcriptomics_analysis/analysis/Jarno## Info: Writing image file mmu04640.pathview.png# Om de afbeelding te weergeven in R:
knitr::include_graphics("mmu04640.pathview.png")
Figuur 30: Hematopoietic cell lineage
Hematopoietic cell lineage: Omdat B-cellen uit hematopoëtische cellen voortkwamen, en de voorloper B-cellen die hiervoor gebruikt zijn nog dichter in de buurt van deze hematopoëtische stamcellen zijn is het volgende pathway ook interessant. Deze pathway maakt ook verschillende dingen zoals de regulatory T-cell en als hier onderexpressie van in zou dit dus invloed hebben op regulatie in de T-cellen. Deze pathway heeft invloed op heel veel essentiele cellen maar omdat nog niet alles bekend is is het moeilijk om te zeggen wat deze over / onder expressies precies doen.
Resultaten & conclusie
Als we kijken naar dit onderzoek hebben we verschillende dingen bekeken, De pathways, De expressie van Ebf1 + Pax5 in jonge en oude muizen en andere genen die hoge expressie hadden. De terug koppeling naar ons originele onderzoek over de verschill van expressie tussen jong en oud van de dus essentiele Ebf1 en Pax5 genen is dat uit onze data is gebleken dat oude muizen lage expressie van deze genen hebben en de jonge muizen juist een hoge expressie, waardoor ik speculeer dat het onderzoek het gelijk had over de A B compartiment (In compartiment A vindt actieve transcriptie plaats en in B niet of in mindere mate.Van Ebf1 is wel bekend dat het met leeftijd van compartiment A naar B gaat, dit komt ook overeen met de verlaagde expressie in de oude muizen (Ebf1) ). Dit zijn compartimenten in de nucleus waar de chromatine heen verplaats waardoor het niet of minder word expressed. In het originele onderzoek was er weinig gespeculeerd over conformatieveranderingen van het DNA dat de onderzochte genen bevat.Wel was experimenteel vastgesteld dat Ebf1 van compartiment A naar B veranderd, op basis daarvan heb ik het bovenstaande vermoeden geconstrueerd.
Er is verder weinig nieuwe informatie gevonden vergeleken met het originele onderzoek ondanks het nieuwe referentie genoom. Veel andere analyses die voor het artikel uitgevoerd waren, waren niet van toepassing voor dit onderzoek omdat het analyses betrok die hier niet gedaan zijn zoals: HI-C, VDJ-seq en HiChIP.